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ABSTRACT 
Land managers frequently need to assess the effects of 
forest cover removal on water yield, water regime, and water 
quality. They may also want to know what happens when forest 
cover is deliberately removed, in a prescribed manner, for 
the purpose of enhancing water supplies. 

One way of approaching these problems is to first 
develop simulation models and then compare model results 
with those obtained in the field. Once the validity of a 
model has been established it may be used to predict the 
impact of a particular forest harvesting pattern on soil 
water and streamflow, or to predict practical watershed 
management prescriptions to improve water yielding 
characteristics of specified river basins. 

Several physically based watershed simulation models 
have been developed for these purposes. Although they 
adequately simulate many hydrological processes, treatment 
of the subsurface flow component is usually less than 
Satisfactory. 

This study is an attempt to rectify the situation 
through development of a physically based, distributed, 
subsurface flow finite element model (SUBFEM), in which 
forest vegetation appears as an integral part. The model 
Simulates the effects of different vegetation patterns on 
soil water distribution and streamflow. Water withdrawal by 
trees is simulated by means of sinks located at appropriate 


near-surface nodes within the finite element mesh. 
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Initial tests on the model, for one- and 
two-dimensional infiltration into originally unsaturated 
soil, produced results that compare favourably with those of 
other investigators. 

Simulations showed that a smaller volume of soil is 
affected by evaporation than by transpiration for a fixed 
evaporation or transpiration demand. Potential gradients are 
greater during evaporation. 

In another set of simulations evapotranspiration from 
Sloping profiles was examined. Trees were assumed to occupy 
different positions on the slope: lower slope only, upper 
Slope only, completely forested, clearcut (no trees), and 
patchcut. 

Quite distinct patterns for both total potential and 
volumetric water content were obtained. They reflected the 
differential drain on soil water by trees and by 
evaporation. 

Time step size, mesh coarseness, and two model 
parameters (initial conditions and saturated hydraulic 
conductivity) were subjected to sensitivity analyses. The 
conclusions drawn from these analyses are outlined below. 

Time step size, under certain conditions, is critical, 
and will determine whether convergence to a solution is 
possible. Mesh coarseness can influence the time at which 
outflow from an originally unsaturated system begins. 

Both the time at which outflow commences, and the 


volume of flow from an originally unsaturated system, are 
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affected by changes in initial conditions. The wetter the 
initial conditions, the earlier outflow begins. SUBFEM is 
very responsive to changes in saturated hydraulic 
conductivity. 

It is felt that the primary usefulness of SUBFEM rests 
in its capability to simulate the processes of infiltration, 
transpiration, and interflow, It can best be used to study 
hydrologic processes on vegetated or partially vegetated 
hillslopes. An example of this kind of application is the 
Simulation of a plot study on Fernow Experimental Forest, 
West Virginia. Data from the field study showed that, for an 
Ppt. Of i. oO. CMAOL rain cover 12 hr, the plot. produced about 
0.4 cm of outflow. The simulation data indicated that the 
entire volume of rain falling on the plot was absorbed by 
the soil. 

Suggestions for improving and applying the model are 


presented. 
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One-dimensional evaporation - total 
potential (a) and water content (b) after 
10,000 sec. 


Finite element discretization for 
two-dimensional horizontal infiltration into 
a an of Yolo Light Clay (After Beven, 
1975). 


Two-dimensional, horizontal infiltration - 
total potential after 6 hr. 


Two-dimensional, horizontal infiltration - 
water content [a) and b)] and total 
potential [c)] fields after 6 hr. 


Two-dimensional, horizontal infiltration - 
advancing wetting front and infiltration 
Capacity. 


Finite element discretization for 
evaporation, transpiration, and 
evapotranspiration from a flat section of 
Yolo Light Clay. 


Evaporation from a flat profile - pressure 
potential (a), total potential (b), and 
water content (c), after 240 hr. 


Transpiration from a flat profile - total 
potential (a) and water content (b) after 
240 hr. 


Evapotranspiration from a flat profile - 
total potential (a) and water content (b) 
after 240 hr. 


Evapotranspiration from a sloping profile - 
total potential (a) and water content (b) 
after 240 hr. 


Evapotranspiration from a sloping profile - 
total potential (a) and water content (b) 
after 480 hr. 


Trees on lower slope - total potential (a) 
and water content (b) after 240 hr. 


Trees on lower slope - total potential (a) 
and water content (b) after 480 hr. 


Trees on upper slope - total potential (a) 
and water content (b) after 240 hr. 
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Trees on upper slope - total potential (a) 
and water content (b) after 480 hr. 


Clearcut simulation - total potential (a) 
and water content (b) after 240 hr. 


Clearcut simulation - total potential (a) 
and water content (b) after 480 hr. 


Patchcut simulation (two cut strips and 
three forested strips) - total potential (a) 
and water content (b) after 240 hr. 


Patchcut simulation (three cut strips and 
two forested strips) - total potential (a) 
and water content (b) after 240 hr. 


Rainfall on a flat profile - total potential 
(a) and water content (b) after 240 hr. 


Rainfall on a sloping profile - total 
potential (a) and water content (b) after 
240 hr. 


Hydrographs for a 10-day rainfall (10 
mm/day), for two mesh sizes, and different 
initial conditions. 


Rainfall on a sloping profile - total 
potential (a) and volumetric water content 
(b) after 240 hr. Obtained using a 217-node, 
360-element mesh. 


Rainfall on a sloping profile with w,;=-10, 
-50, and -100 cm - total potential (a) and 
volumetric water content (b) after 240 hr. 
Fine mesh simulation. 


Recession curves - obtained when rainfall 
rate is changed from 10 to 2 mm/day, and 
time step sizes are 0.1 and 0.01 hr. 


Recession curves - obtained when rainfall 
rate is changed from 10 to 2 mm/day, and 
time step sizes are 6, 3, and 1 hr. 


Recession curves - obtained when rainfall 
rate is changed from 10 to 0 mm/day, and 
time step sizes are 0.1 and 0.01 hr. 


Recession curves - obtained when rainfall 


rate is changed from 10 to 0 mm/day, and 
time step sizes are 6, 3, and 1 hr. 
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Total potential at end of first time step 
following cessation of 10 cm rainfall. 


Hydrographs 
sizes using 


Hydrographs 
coarseness, 


Hydrographs 
coarseness, 


Hydrographs 


obtained for different time step 
a 64-node, 90-element mesh. 


produced by meshes of different 
and a time step size of 24 hr. 


produced by meshes of different 
and a time step size of 3 hr. 


obtained using a 114-node mesh 


and different initial conditions. 


Variation of discharge over time for 
different values of saturated hydraulic 
conductivity (K,). 


Water content in Yolo Light Clay after 10 
days of rainfall, for two values of 
Saturated hydraulic conductivity (K,). 
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I. INTRODUCTION 


The eastern slopes of the Rocky Mountains in Alberta 
form the headwaters of the Athabasca and Saskatchewan River 
systems. The Saskatchewan River, which drains mountains, 
foothills and plains, is the primary source of water for the 
Canadian Prairies, and is therefore extremely important to 
the economy of that area. A far greater proportion of the 
Saskatchewan River's total annual flow is derived from the 
heavily forested mountains and foothills than from the 
plains. Consequently, water production is considered by some 
to be the prime use of the East Slopes region. In the past, 
management of the Bast Slopes watersheds has been 
essentially protective in nature. However, as population 
pressures increase, the demands for good quality water also 
increase and ways must be found to meet these demands. One 
way of increasing the water supply is to remove, ina 
prescribed manner, some of the forest cover from the high 
water-~yielding portions of the watershed. 

Industrial activity such as forestry, coal, oil, and 
gas exploration is increasing on the East Slopes as the 
search for, and the exploitation of, additional natural 
resources intensifies. It is important that land managers be 
able to assess the effects of such activity on water yield, 
water regime, and water quality since there is convincing 
evidence that removing forest cover for wood production or 


other industrial uses can result in an increase in stream 
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discharge from the cutover areas. A reduction in water 
quality can also result if due consideration is not given to 
proper design of roads and water course crossings. 

One way of approaching the kind of problems outlined 
above is to first develop watershed simulation models and 
then compare model results with those obtained in the field. 
Once the validity of a model has been established it may be 
used to predict the impact of a particular forest harvesting 
pattern on streamflow. Alternatively, it can be used to 
predict practical watershed management prescriptions to 
improve water yielding characteristics of specified river 
basins. 

The type of conceptual model envisaged here is one in 
which the individual components or hydrologic processes are 
related to various forest stand parameters such as species 
composition, age, stand density (stems per hectare), basal 
area (m? per hectare), and areal distribution. When all 
trees are removed from an area, Such as occurs during 
clearcutting, all the stand parameters for that area become 
equal to Ss 

For many hydrologic processes, models based on these 
kind of relationships and which produce satisfactory results 
are available. They are well documented in the literature. 
For example, several models have been developed which 
Simulate the effects of the forest and forest harvesting on 
Snow accumulation, redistribution, and melt. Other processes 


as related to forest stand parameters have not been so well 
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defined. Subsurface flow is perhaps the most complex of 
these. 

Ideally, the preferred forest vegetation manipulation 
watershed model is one which is completely physically based 
and which is capable of predicting with a fair degree of 
accuracy: a) the cutting pattern and amount of forest to be 
removed in order to obtain a hydrograph with specified 
characteristics, b) how a variety of land uses will change 
the hydrographs of affected streams, and c) the changes in 
water quality produced by forestry and other industrial 
operations on the watershed. The model must be simple enough 
that it can be applied by forest, land, and water managers 
to areas for which data is not extensive. A watershed model 
possessing all these characteristics does not exist. 

Although many watershed models are available, few are 
designed to solve forest watershed problems, and fewer still 
contain a physically based subsurface flow component. Often, 
this component is simulated, if at all, in a somewhat 
arbitrary manner. 

The purpose of this project is to progress toward the 
idealized model described above, by developing a physically 
based subsurface flow medel which provides for simulation of 
forest stands. The subsurface model will: a) identify the 
main flow paths taken by water through the porous media to 
the stream channel after it passes through the air-soil 
interface, and b) quantify the outflow from the seepage face 


at the stream channel, so that streamflow response to 
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watershed treatment and precipitation can be determined. 

It is felt that this purpose can be best achieved by 
pursuing the following objectives: 

1. to develop from Darcy's Law and the continuity equation 
a mathematical model of two-dimensional transient 
unsaturated and saturated flow through porous media, 
applicable under natural conditions; 

2. to incorporate the mathematical model as part of a 
physically based synthesis of the hydrologic cycle; 

3. to use the synthesized model to simulate behaviour of 
hillslopes following treatment, and to check simulation 


results against field data from treated hillslopes. 
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II. LITERATURE REVIEW 


A. History Of Forest Hydrology And Watershed Management 
Research 

The interrelationships between forests and water have 
long been recognized. As early as 1215, Louis VI of France 
issued "The decree of water and forests." In 1342, a 
community in Switzerland reserved the first protection 
forest as a safeguard against avalanches (Kittredge, 1948). 
Since that time many countries have taken similar and 
additional measures such as reforestation and construction 
of check dams to curtail erosion or to prevent damage to 
life and property by torrents and avalanches. However, most 
of the significant eet ent tite and legal developments related 
to forest-water interactions have occurred over the past 150 
years. 

During this period the United States emerged as a 
leader in one field of forest-water relations research and 
many reports were published on the effects, deleterious or 
otherwise, of forest cover removal on climate, soil and 
water. In 1863 for example, G.P. Marsh wrote a book entitled 
"Man and Nature", which was later revised and republished 
several times under the title "The Barth as Modified by 
Human Action". This book created a national awareness toward 
watershed management. Several years later Kinney (1900) 
discussed the problems associated with water, forest fires, 


and exploitation of forest lands in Southern California. 
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Such reports had a significant impact on U.S. land policy 
formulation and on ensuing legislation. 

It is evident from the wording of legislation 
authorizing the reservation of forested public lands that 
.the water-controlling function of the forest was to be 
regarded as equally important as its capacity to supply 


timber. Thus, one clause in the Act of June 4, 1897 reads 


",.-no public reservation shall be established 
except to improve and protect the forest within the 
reservation or for the purpose of securing favorable 
conditions of water flow and to furnish a continuous 
supply of timber." (Kittredge, 1948). 


Some of the more important pieces of U.S. Federal 
legislation together with the names of the key people 
involved have been summarized concisely by Hewlett and 


Nutter (1969) 


"A period of evaluation and alarming descriptions 
(Marsh's' book, 1863) led to the propaganda period 
of forest influences (1877-1910). Forest influences 
played a large role in forestry, conservation and 
public land policies. Fernow, Hough, Pinchot, Roth) 
and Roosevelt helped to formulate policies which 
produced many acts of Congress related to forests 
and water: Weeks Law (1911); Forest Experiment 
Stations authorized by McSweeny-McNary Act (1928); 
New Deal of the 30's, including the CCC Program; The 
Flood Control Act (1936); Soil Conservation Service 
(1936). Public Law 566 (Small Watershed Act of 1954) 
forced involvement of local people in watershed 
management. The Senate Select Committee (1958-61) 
found that training, planning, and pollution were 
the greatest water problems. Since then a flood of 
acts aimed at water resources planning, 
conservation, and development have been passed. 
Several important ones were the Water Resources 

' Marsh, G.P: 1863. Man and Nature. Scribner and Sons, New 

York, NY. (Later titled Earth as Modified by Human Action). | 
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Research Act (1965), Water Pollution Control Act 
: Ge and the Water Resources Planning Act 
1966)." 


In Canada, forest-water relations research has not 
enjoyed the recognition accorded this field in the United 
States. Nevertheless, its importance was recognized as early 
as 1910 when legislation was passed to establish the Rocky 
Mountains Forest Reserve in the Saskatchewan River 
headwaters of Alberta as a watershed protection zone 
(Alberta Energy and Natural Resources, 1979). 

Several collective terms have been used to describe 
forest-water relations. In older literature, the term 
"forest influences" is used. Forest influences have been 
defined to include all effects resulting from the presence 
of forest or brush upon climate, soil water, runoff, 
Streamflow, floods, erosion, and soil productivity 


(Kittredge, 1948). It has been replaced by "forest 


hydrology” which is 


"a branch of hydrology that deals with the effects 
of forests and associated wildland vegetation on the 
water cycle, including the effect on erosion, water 
quality, and the microclimate" (Hewlett and Nutter, 


969°)’. 


This term has been used interchangeably with "watershed 


management" which has been defined as 


"Application of business methods and technical 
principles to the handling of all the renewable 
resources in a watershed to assure maximum supplies 
of useable water, desirable waterflow, prevention 
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and control of erosion, and the reduction of flood 
and sediment damages.” (Society of American 
Foresters, 1958). 


An alternative definition of watershed management is: "The 
management of land for optimal production of water, with due 
attention being given to soil stability and to the other 
resources of the land." (Jeffrey, 1969). 

Generally speaking, "forest hydrology" refers to the 
Scientific-technical basis for forest watershed management, 
-and "watershed management" refers to the management of land 
for water production itself (Jeffrey, 1969). Other terms 
such as "wildland hydrology" and "land use hydrology" are 
also used, but in the forest environment context they are 
Synonymous with "forest hydrology". 

Forest hydrology then can be considered as the study of 
hydrologic processes such as precipitation, interception, 
overland flow, infiltration, evapotranspiration, groundwater 
flow, streamflow, erosion, sedimentation and water quality 
in a forest setting. Its subject matter also includes study 
of the side effects of forestry operations (timber 
harvesting, regeneration, tree planting, vegetation type 
conversion, and the application of fertilizers and 
pesticides), fire and grazing on water supply, floods, 
erosion, and water quality. 

The scope of watershed management includes management 
of the forest land for one or more of the following specific 


objectives: optimum water yield, maximum water yield, flood 
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reduction, and maintaining water quality? by minimizing 
erosion and stream enumen tan car Typically such management 
is conducted in high- to medium-water yielding source or 
headwaters areas of river basins. Generally such areas are 
located in mountainous or hilly terrain, and are subject to 
greater precipitation than the plains below. 

Forest hydrology and forest watershed management 
research have been pursued vigorously during the twentieth 
century, particularly in the United States. An excellent 
Summary of American experience and results is provided by 
Anderson et al. (1976). Several important textbooks have 
been published both in North America (Kittredge, 1948; 
Colman, 1953; Hewlett and Nutter, 1969; and Satterlund, 
1972) and in Europe (Geiger, 1950; Molchanov, 1960; and 
Rakhmanov, 1962). Several manuals (Food and Agriculture 
Sncanizati oncofuthéeUniteduNations-;ni97b6a, o19¢6byeng77; U.S. 
Forest Service, 1974) were made available, to assist the 
practicing forest manager. Two important international 
symposia on forests and water were also held (Sopper and 
Lull, 1967; Talat and Dunford, 1970). Several other symposia 
were organized at the local level (American Society of 
Agricultural Engineers and American Society of Civil 
Engineers, 1970; Csallany et al., 1972; Krygier and Hall, 
1971: Society of American Foresters and Oregon State 


2 Maintaining water quality through erosion and sediment 
control is an extremely important objective in watershed 
Management and is the primary function of many protection 
watersheds. However, the subject matter 1s beyond the scope 
of this paper and will not be discussed in detail. 
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University, 1963, 1966). Another international symposium - 
on the results of research on representative and 
experimental basins - was held in Wellington, New Zealand 
(International Association of Scientific Hydrology and 
UNESCO, 1970). Although it did not deai exclusively with 
forests/water problems, the subject matter of the symposium 
was highly relevant. 

Tables 1, 2 and 3 are taken from Anderson et al. (1976) 
and from Hibbert (1967) and summarize the results from much 
of the research referred to above. Table 1 shows the 
increases in water yield following forest cutting in various 
parts of the United States and in other countries. 

The Emmenthal watershed study begun in 1900 predates 
the studies referred to in Table 1 and involved two small 
watersheds in the Emmenthal Mountains of Switzerland. One 
watershed was completely forested and the other mostly 
pastureland. The objective was to determine the influence of 
the forest on the water balance. Although differences (246 
mm) in streamflow between the two watersheds were detected, 
ieawas mot. possible io attribute these, diffterencess tos the 
influence of forest cover alone (Colman, 1953). However, the 
Study was the first attempt to solve a forest hydrology 
problem using entire watersheds and scientific methods. 

The Wagon Wheel Gap experiment in Colorado (Table 1) 
was begun in 1910 and is significant because it was the 
first time that the control (paired) watershed method was 


used. This method requires measurement of streamflow from a 
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control sub-watershed and from an adjacent to-be-treated 
sub-watershed. Flow data are collected for a calibration 
period of, preferably, 10 or more years, and a statistical 
relationship obtained between the two sets of flow data. A 
treatment such as commercial logging and associated road 
construction is imposed on one sub-watershed. This marks the 
end of the calibration phase and the beginning of the 
post-treatment period. Flow data are obtained from both 
sub-watersheds for some minimum period following treatment, 
and a statistical relationship obtained. If the relationship 
1s Significantly different from that determined for the 
calibration period data, then this is regarded as evidence 
that the treatment has affected streamflow. This method was 
to be the basis for much of the forest hydrology research 
that followed. For Wagon Wheel Gap, the investigators showed 
conclusively (Fig. 1) that removing 100 percent of the aspen 
and Douglas fir resulted in a measurable increase in water 
yield. 

Data (Table 1) from two different watersheds will be 
used to illustrate the great variability in results due to 
geographic location, climatic and other factors. In 1940, 
one subbasin of the Coweeta watershed in North Carolina, 16 
ha in extent and supporting mixed hardwoods, was 100 percent 
clearcut. Located in the humid mountain region of the United 
States, it has a mean annual precipitation (all rain) of 
1829 mm and a mean annual streamflow of 787 mm. During the 


first five years following cutting, water yield increases 
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were recorded as 366, 277, 277, 249, and 201 mm. The first 
value represents a 66 percent and the last a 31 percent 
increase in water yield. Evidently, the effect diminishes as 
regrowth occurs. The experiment was repeated in 1962 and the 
results for the first year following treatment were 
identical to the corresponding year of the earlier 
experiment (Hibbert, 1967). 

The Fraser watershed study in the high-elevation zone 
of the Colorado Rocky Mountains consists of Fool Creek and 
East St. Louis Creek where the latter serves as control. 
Mean precipitation for the area (Table 1) is 762 mm, 75 
percent of which occurs as snow. The mean annual streamflow 
is 279 mm. In 1954, 40 percent of the lodgepole pine, 
Spruce-fir forest on Fool Creek watershed, which is 289 ha 
in area, was clearcut. Alternate clearcut strips 20, 40, 60, 
and 121 m wide were extended up and down slope and bounded 
at the ends by contour roads. This treatment produced 
increases in water yield of 84, 132, 94, 117, and 137 mm 
during each of the next five years. These figures are 
considerably lower than those from Coweeta, but the highest 
values for both Coweeta and Fool Creek are about 18 to 20 
percent of mean annual precipitation. In the Colorado study 
the smallest increase (32 percent) in water yield occurred 
during the first year following cutting and the largest (71 
percent) during the fifth year. The effect is not diminished 
over a short time period as is the case in the Coweeta 


experiment. In fact, increases in flow of 51 to 102 mm from 
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Fool Creek have been sustained for 17 years (Fig. 2), and 
are expected to persist for considerably longer (Anderson et 
al. 1976). The main factor at Fraser is the considerable 
length of time required for coniferous subalpine forests to 
grow to maturity. In contrast, at Coweeta where the increase 
in water yield diminishes rapidly with time, the climate is 
more favorable to rapid regrowth of forest cover. 

The mechanism by which increased water yields are 
generated differ between the two locations. In Coweeta extra 
water becomes available through reduction in transpiration. 
In Colorado, the primary cause for increased water yields 
from subalpine forests is the combined action of snow 
redistribution and reduced evapotranspiration (Leaf, 1975). 

Table 2 summarizes results from studies investigating 
the effects of reforestation and afforestation (establishing 
a forest where no forest existed before). In nearly every 
case there is a reduction in water yield. The decreases 
range from 0 to 277 mm and appear to be effective for as 
long as the new forest cover exists. Factors contributing to 
water yield reduction may include: increased interception 
and infiltration rates, diminished overland flow, increased 
detention storage, and increased transpiration. 
Reforestation may also lead to reduction in peak flows 
(Table 3) from a watershed. The reduction varies with the 
type of cover that existed before conversion, the proportion 


of the watershed planted, and with the season of the year. 
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The effects of forest cover removal on peak flows are 
still unclear. Anderson et al. (1976) give details of many 
studies where this factor has been investigated and they are 
repeated here. Two examples are taken from Coweeta. The 
first 1s a 13 ha, 100 percent clearcut watershed where no 
increase in maximum peak discharge was detected. The second, 
a 44 ha watershed was also clearcut 100 percent, and the 
peak flows increased by 9 percent. 

A study of commercial logging on watersheds in West 
Virginia showed that the effects on storm peaks depended on 
the season. Eighty-six percent of the total wood volume was 
removed from a hardwood-forested watershed. Instantaneous 
peaks during the growing season increased by an average of 
21 percent; in the dormant season they were apparently 
reduced by 4 percent. In Japan, clearcutting a 2.4 ha 
watershed increased peak runoff from heavy rains by about 20 
percent. 

Both forest and ground cover were removed from the 
Hubbard Brook Experimental Forest in New Hampshire. Summer 
peak flows increased considerably. For the six highest peak 
flows during June through September, 1966 through 1969 the 
peaks for a 16 ha denuded watershed averaged double the 
expected amount, with eeeaees for the individual events 
ranging from -19 to +250 percent. 

Clearcutting can either increase or decrease peak flow 
rates when snowmelt is involved. In another experiment at 


Hubbard Brook, clearcutting a small watershed increased peak 
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flows early in the snowmelt season by as much as 35 percent 
and decreased them as much as 66 en eene later in the 
season. 

At the H.J. Andrews Experimental Forest in central 
Oregon during the 1965-69 period, the peaks for 84 percent 
of the storms were greater in the clearcut watershed than in 
the uncut control. Other studies show that clearcutting 
watersheds in Oregon may produce an increase in peaks from 
fall storms of 90 percent and from winter storms of 28 
percent. (Anderson et al. 1976) 

At Fraser, Colorado the treatment referred to earlier 
produced no significant change in the average peak flows 
(Leaf, 1975), although peak discharges from snowmelt were 
greater than predicted, increasing by 50 and 45 percent 
during the first and third years following treatment, 
respectively. The peak for the second year was 23 percent 
less than had been expected. (Anderson et al. 1976) 

Given the results from the investigations described 
above together with results published elsewhere, it is 
possible to illustrate most of the facets of forest 
hydrology and watershed management in the form of a 


reversible quasi-equilibrium equation: 
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Impact e Disaster - fire, insect, or disease 
2 Harvest 
Manipulative 3 Watershed Management 


—_—_—-_———————————— WATER 


Processes 
(Black box) 
FOREST a 


4 Reforestation - planting, or natural 
5 Afforestation 


The events shown above the equation result in reduced forest 


cover and increased water yield or a shift in equilibrium to 


the right. For events listed below the equation, the reverse 


aS tinue. 


It will be observed that the watersheds referred to so 


far are all quite small - only one exceeds 500 ha in area. 


Therefore, the results to be summarized are applicable 


primarily to small watersheds. It appears that removing 


forest cover affects hydrological processes in the following 


ways: 


decreases interception (no forest canopy); 

reduces transpiration (no trees); 

increases evaporation (reduced ground protection); 
reduces soil water deficits (decreased draft on soil 
water); 

increases water yield, the absolute increase being 
greater in humid than in dryer regions; 


increases peak flows; 


ae 
Fa 


rO2 7 U 
’ ¢ ‘ 
‘ c or 
vt 
4 
** 4 
: 1 
105 Goegs 242 
= i e* 
A Ahan =~ 


3 7s ~ 4 
tes 
=) « bu 
eerbeh.) 
e4nfasds 


iy 


[a 


as 


4 

- 5 ee a 
e PI 

a a ~ 


- 
= 


nf Sekto) teal 
88 One® a s 


24 


7. changes the aerodynamic effect of the forest (wind flow 
patterns and velocities are altered which has important 
implications for snow redistribution); 

8. in areas where snow forms a significant part of the 
precipitation, increases snow accumulation in the open 
areas and reduces it in the adjacent forest. (Snow in 
the open is also subject to increased sublimation. Under 
certain sheltered conditions, snow tends to remain 
longer in the openings than in the adjacent forest); 

9. water yield increases are uSually greatest during the 
first few years following forest cover removal. The 
effect is reduced as the area becomes restocked with 
trees. In humid areas where plants grow rapidly, the 
increase in water yield may be significantly reduced 
within a few years after forest clearing. In areas where 
snow is an important factor and where plant growth is 


slow, the effects may remain for more than 30 years. 


When regrowth is established, or when reforestation or 
afforestation are effected, the foregoing processes are 
reversed. This is a generalization, but there are usually 
increases in interception, transpiration, infiltration rate, 
and detention storage, as well as reductions in water yields 
and peak flows. As the forest canopy becomes more uniform in 


height the snow depth under the canopy also tends to become 


more uniform. 
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Another important finding from many studies is that 
Hortonian overland flow is not a significant process on 
forested watersheds except on and in the vicinity of logging 
roads and exposed soils. Usually infiltration capacity of 
soils under forest exceeds the rainfall intensity by a 
comfortable margin. 

In spite of the progress made in the fields of forest 
hydrology and watershed management, many problems remain. 
Most important perhaps is the insufficient attention paid to 
the necessity of extrapolating results from small watersheds 
to large ones for land management purposes. The problem is 
compounded by contradictions between results reported for 
small watersheds and those obtained for large river basins. 
Russian investigators (Rakhmanov, 1962, 1970a, 1970b; and 
Bochkov, 1970) used statistical procedures to show that a 
positive relation exists between percent forest cover and 
Streamflow for large river basins. In the case of the Upper 
Volga Basin study (Rakhmanov, 1970b) this conclusion was 
based on data (Fig. 3) from 53 basins ranging in size from 
226° tor 13,528 ‘km*. However, the positive relationship 
discovered by Rakhmanov and Bochkov exists because more 
trees grow where water and climate are favorable; it does 
not imply that water yields are reduced when forest cover is 
decreased. 

There are two instances where results from large basins 
have supported conclusions drawn from small watershed 


Studies. The first concerns destruction of a forest by 
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insects (Love, 1955), and the second involves long term 
recovery of forest vegetation in the Sacandaga River 
drainage, New York (Hibbert, 1967). An outbreak of Engelmann 
Spruce beetle on the White River Plateau of western Colorado 
killed Engelmann spruce and lodgepole pine on 30 percent of 
a 1974 km? drainage of the White River. The increase in 
average annual streamflow of 58 mm during the years 
following the outbreak was attributed to reduced 
interception of snow and reduced transpiration by 
beetle-killed trees (Love, 1955). Results from the New York 
Study (Hibbert, 1967) showed that increasing the forest 
basal area from 17 to 28 m?/ha on the 1272 km? Sacandaga 
River watershed over 38 years caused a reduction in water 
yield of 196 mm. 

Although a considerable amount of information is 
available on the behavior of watersheds, it is a remarkable 
fact that there are no examples of watersheds, large or 
small, being managed for increased water production. Forest 
watershed management today is almost exclusively protective 
in nature, directed primarily to preventing erosion and 
maintaining water quality. However, some estimates are 
available on the potential for increasing water yield in the 
United States (Table 4). 

From the foregoing discussion it is evident that 
considerable research still remains to be done and should be 
directed to solving the following problems: 


1. develop extrapolation procedures so that results from 
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small experimental watersheds can be applied to 

a. Similar basins in the same region, 

b. basins with different characteristics, 

c. large river basins; 

in the absence of Hortonian flow, obtain a better 
understanding of subsurface flow as influenced by 
geology, vegetation, landform, and other important 
factors; 

develop improved methods of routing flow through the 


Subsurface portion of the watershed. 


Watershed experiments are costly and time-consuming 


undertakings. Because of this, and the availablility of 


powerful computers, researchers are turning increasingly to 


watershed simulation models for answers to their watershed 


problems. Nevertheless, gauged watersheds are still 


necessary to provide field data against which model outputs 


can be checked. This thesis is an attempt to find solutions 


to the subsurface flow problems described in 2 and 3 above, 


in the particular setting of a forested foothill or mountain 


watershed, using computer simulation models. 


B. Development Of Forest Hydrology And Watershed Management 


In Alberta 


The Saskatchewan River system which rises on the east 


slopes of the Rocky Mountains of Alberta is an important 


Source of water for the prairie region of Alberta, 
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Saskatchewan and Manitoba. About 87 percent of the drainage 
area is located on the plains, the remainder lies within the 
foothills and mountains of Alberta. However, this smaller 
portion (13 percent) contains the upper headwaters of the 
Saskatchewan River system and contributes 87 percent of the 
total annual flow (Hanson, n.d.). 

The importance of these eee headwaters was 
recognized as early as 1910 when the Dominion Government 
established the Rocky Mountains Forest Reserve, an area of 
about 23150 km*. The Forest Reserve contains most of the 
Saskatchewan River headwaters area not included in the 
National Parks (Fig. 4). 

Later, discussions between the Canadian Government and 
the Province of Alberta led to the formation in 1947 of a 
joint federal-provincial agency, the Eastern Rockies Forest 
Conservation Board. Its task was two-fold: 

1. to operate a program designed to protect forests from 
fire, insects, disease and other damage; 

2. to conserve, develop and manage forests to obtain the 
optimum flow of water (optimum quality, quantity and 
timing) in the Saskatchewan River and its tributaries 


(Alberta Energy and Natural Resources, 1979). 


After a decade of operation, the Board faced several 
problems: its ability to function effectively was hampered 
by lack of basic data from research, local water shortages 


were becoming more common, and there were conflicts in 


Sechrest sd? To) JAsoIeg ‘VE mud & sori Taw antag 
a4 niddiw Sebhl iShlemet ark ere ais noha 


on 
a 
«sf tewe oidts ,asvegoM .4?2sela ie biainiaeeaall 


; . ‘on 
ons Ys atesaW¥tast wsd0gy sit @Arszuon (Iago seq er) f 
bne sstéys Saeke 
bn , ose Wee © 


ow DESEJIEVRSE A Bega. ic. aweott. 3I£ ones tog 


frien ts 708 im tf etioty gaawy 0 4 @acé <i :668 ae 
q 


+) Bets AB YeNUIgeRs Seeusl joi eorio Wane wes 
gf= 20 t26m anise 4¥tegen sze 2% 4) . eR 


e439 nit BaeBiiegen: 40 gence = 44e4¥Oe0n THRSF Te 


2 . 29) 


m2 Peis ret PISS 8 ores -{- & 
» As ant 
thbsoaét). Taiscn Ine e&ese2y ec somith ae « 

re 


atu nisidh.icg-e Seats? oganeg fia6 qola-et ov TRsHOD r 
al - 7 : . 


r : , i 
Bie Yori insu yesiler, Misisgs) serev 36 “oi: corn i 


garrewatwa. 2et Bue sevict cewad>tadase oft Wh tt 


7 f 
vA > she t _ 
ry . (2Eer A fesu7e% bas yideal” 
_ 


of _ fy - Ea ; 
Leisver Dea ® Jas ane ,WetFexegc to sba oo 0 } até 
i, ae is oe i ° 


ae negssaut 93. re if ide & 


aw 


a siltieeess aor4 ash: 


id, 


a 


- ALBERTA~ 


S) 


50 75 100 km 


25 fe} fs 50 75 miles 


cm oo 


"=; CANADA 


UNIT, STATES 


Figure 4. Map showing Rocky Mountains Forest Reserve (Adapted from 
maps published by the Alberta Government). 
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demand for available water supplies downstream (Redmond, 

1964). 

In 1959, senior federal and provincial authorities met 
to examine these problems. They agreed that a special 
Committee - the Technical Advisory Committee - be formed to 
define research and management problems related to watershed 
management on the east slope of the Rocky Mountains. This 
Committee would carry out a research program which would 
provide: 

1. a better understanding of hydrologic processes both 
locally and in the entire area; 

2. evaluation of newly developed and old methods of 
protecting and maintaining watersheds while forests and 
other resources are being utilized; 

3. means of restoring damaged watersheds; 

4, methods of improving streamflow (quantity, quality, and 


timing) from the watershed. 


This program came to be known as the East Slopes 
(Alberta) Watershed Research Program. It was also agreed 
that most of the research should be carried out within the 
Rocky Mountains Forest Reserve, and that small 
representative watersheds be selected for intensive study 


(Redmond, 1964). 
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Jeffrey*® was perhaps the first to attempt a problem analysis 
and to suggest methods by which stated research program 
objectives might be attained. Subsequently, three 
representative basins were selected for study. Marmot Creek 
Basin (940 ha in area) was established to study the 
hydrology of spruce-fir forests. Deer Creek Basin (544 ha) 
was selected as representative of lodgepole pine forest, and 
Streeter Basin (598 ha) was chosen for research into montane 
aspen forest and associated grasslands. These three 
vegetation types comprise most of the Forest Reserve area 
which lends itself to any type of vegetation manipulation 
for the purpose of altering water yield and regime. 

Deer Creek Basin was later withdrawn from the Watershed 
Research Program because the commercial harvest objectives 
set for the lodgepole pine basin had been realized elsewhere 
on the East Slopes. It was established that harvesting 
lodgepole pine and spruce-fir does increase streamflow 
(Swanson, 1977; Swanson and Hillman, 1977). 

During the decade 1960 to 1970 meteorological, 
groundwater, hydrometric and other instrumentation networks 
for the three basins were completed (Jeffrey, 1964, i965; 
East Slopes (Alberta) Watershed Research Program, 1966). In 
1969, when the Program was expanded to include the entire 
province of Alberta, its name was changed to Alberta 
Watershed Research Program (AWRP) (Swanson, 1977). The AWRP 
Jeffrey, W. W. 1961. Prerequisites and Priorities for 
watershed research in the Eastern Rockies, Alberta - a 


tentative initial appraisal. Can. Dep. For., For. Res. Br., 
Calgary, Alberta. Unpubl. rep. 
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has become the primary vehicle for conducting and overseeing 
forest watershed research in Alberta. It replaces the 
Eastern Rockies Forest Conservation Board which ceased to 
exist in 1973. 

The expansion introduced other watershed experiments, 
located outside the Forest Reserve, into the program. 
Tri-Creeks Basin for example is being monitored to determine 
the impact of a pulp timber harvesting operation on 
lodgepole pine-foothills hydrology with special attention 
being devoted to maintaining water quality and preserving 
fish habitat. Spring Creek project was established to 
examine the hydrologic effects of converting aspen-Sspruce 
forest land for agricultural use. 

The principal objectives for each basin are similar to 
the objectives of the AWRP as a whole: to establish the 
hydrology of the basin, to assess the effect of forest cover 
removal on water quality, yield and regime, and/or to 
develop vegetation manipulating schemes which will produce 
favorable changes in ee quality, yield or regime. The 
common procedure adopted to attain these objectives is the 
paired watershed approach described earlier. 

The program's first treatment, a "commercial" forest 
harvest, was carried out on Cabin Creek Basin, a subbasin of 
Marmot Creek. Another subbasin, Middle Creek Basin, which 
adjoins Cabin Creek subbasin, serves as control. Logging 
haul roads were constructed in Cabin Creek Basin during 


1971-72. The commercial harvest itself was not Carried out 
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until 1974. Data for evaluating the effects of road 
construction alone on sediment production in Cabin Creek 
were obtained during the intervening period. The commercial 
harvest on Cabin Creek Basin was completed in October, 1974. 
Six blocks representing 50 percent of the forested area in 
the subbasin were clearcut. 

In 1976, the west subbasin of Streeter Basin was 
treated to improve range conditions for cattle and wildlife 
while maintaining or prolonging streamflow during the summer 
months. The hydrologic objectives were to be met by altering 
the vegetation so that a) snow would be concentrated in 
areas where sublimation and melt would be minimized, and b) 
infiltration of snowmelt would be enhanced by reducing snow 
accumulation in discharge areas and increasing i eS eo 
recharge areas. 

For this purpose a total area of about 55 hectares was 
treated. Patches 1 1/2 to 2 tree heights in width and 60 to 
180 m in length were cut in the aspen stands. In shrub 
stands, strips of 3 to 5 m width were cleared (Golding, 
177). 

A second treatment for Marmot Creek Basin commenced in 
1977. In this case the objective is to alter the water 
regime or, more specifically, to delay the time of peak 
runoff and to prolong recession flow resulting from 
Snowmelt. Middle Creek subbasin again serves as control. 
Treatment was carried out on Twin Creek subbasin and 


consisted of cutting 40 percent of the forested area in 2500 
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circular openings. The size of openings ranged from about 
3/4 to 1 1/2 times the height of the adjacent trees 
(Golding, 1977). Treatment was completed in late 1979. 

In August 1978, the Oldman River Basin Study Management 
Committee recommended that research areas be established in 
the Oldman River Basin, Southern Alberta to provide 
convincing data on whether vegetative manipulation would 
produce beneficial increases in water yield (Environmental 
Council of Alberta, 1979). Unlike the experiments described 
previously, this study would be conducted over a large area 
(> 260 km?). The recommendation calls for tests of three 
different forest land management practices: 

1. no vegetation manipulation; 

2. harvesting using good commercial harvesting techniques 
as approved by the Alberta Forest Service; 

3. vegetation manipulation for the sole purpose of 
increasing water yield. These practices would be 


monitored on three representative sub-basins. 


The recommendation is an outcome of a report, prepared 
for the Oldman River Basin Study Management Committee 
(Northwest Hydraulics Consultants Ltd., Hyat Resource 
Services Ltd., Edmonton, Alberta. 1977) which examined the 
possibility of increasing water yields through watershed 
Management in the upper Oldman River Basin. This report is 
one of several commissioned by the Committee, in response to 


water shortages in Irrigation Districts of southern Alberta, 
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to investigate alternatives for increasing water yields in 


the Oldman River Basin. 


C. The Watershed System 

Forest watershed management and forest hydrology 
Studies are usually conducted on small watersheds that 
provide flow during most of the year and are less than 25 
km? in extent. The main reason for selecting small 
watersheds is that travel time and channel storage during 
storm periods are not the major factors in determining 
volume of direct runoff. If they were, then the land-use 
effects being measured would be masked by the channel 
Storage effect. 

A watershed system consists of the area contained 
within topographic divides, the precipitation falling on the 
area, water stored or detained by the watershed, and the 
discharge from its exit point. Vertically, the system 
extends high enough above the vegetation canopy to include 
both total precipitation and the aerodynamic effects of the 
canopy. It extends low enough to include all porous media 
and subsurface water influencing or contributing to channel 
flow within the watershed. Usually, it 1S assumed that the 
groundwater or phreatic divide coincides with the 
topographic divide, and that water enters the watershed as 
precipitation and leaves it as channel flow or is returned 
to the atmosphere. The second assumption implies that there 


is no subsurface inflow and outflow i.e. the watershed does 
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not leak. Such assumptions must be verified by appropriate 
subsurface investigations. 

The hydrological processes governing water movement 
through a watershed are well documented. Many processes such 
as infiltration, flow through porous media, and open channel 
flow can be expressed in the form of one or more equations, 
the derivation of which is based on sound physical 
principles. Most hydrologists accept these relations when 
the processes are considered in isolation, but often 
disagree on how these processes should be integrated to 
describe the watershed system. A prime example of such a 
Subject of disagreement is the movement of water from the 
land surface to the stream channel. 

Hydrologists have long recognized that water moving 
from the ground surface to the stream channel follows three 
routes of travel: overland flow, interflow, and groundwater 
flow. Overland flow or surface runoff is that water which 
travels over the ground surface to the channel. Interflow is 
water which infiltrates the soil surface and moves laterally 
through the upper soil layers. Groundwater flow, also called 
base flow or dry-weather flow, occurs when the water table 
intersects the stream channel. Usually, the total flow is 
divided into two parts: storm, or direct, runoff and base 
flow. Direct runoff consists of overland flow and interflow, 
and baseflow is groundwater. This distinction is based on 


time of arrival in the stream channel rather than on the 


path followed (Linsley et al., 1958). The relative 
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importance of overland flow, interflow and groundwater flow 
has been the subject of much debate among hydrologists. 

The classical concepts of infiltration and overland 
flow were first propounded by Horton (1933, 1936, 1945) 
nearly 50 years ago. He developed a negative-exponential 
decay function which determined the changes in infiltration 
Capacity aS a rainstorm progressed. Simply put, Horton's 
theory states that if the rainfall intensity exceeds the 
infiltration capacity, then the excess rainfall becomes 
overland flow which as the sole determinant in the stream 
hydrograph peak. Overland flow is assumed to occur over all 
parts of the watershed as a thin film or sheet. 

These concepts were not seriously questioned until the 
1960's when several investigators presented new insights 
into process and watershed behaviour. Hewlett (1961) 
observed no overland flow on a forested watershed at 
Coweeta, North Carolina. Further investigation also showed 
an absence of any large groundwater system to provide 
nonstorm streamflow. Using measurements of saturated and 
unsaturated flow through a steeply-inclined trough, Hewlett 
(1961) showed that unsaturated soil can serve as main 
storage or source for base flow. He concluded that: 
"up-slope rain charges the soil mantle in preparation for 
Succeeding days and weeks of base flow, whereas downslope 
rain and channel interception will furnish most of the 
Streamflow". Thus the areas immediately adjacent to small 


streams provide, through groundwater drainage, a 
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disproportionate percentage of stormflow relative to the 
total area of the watershed. This implies that, contrary to 
Horton's theory, only a small percentage of a watershed area 
contributes storm runoff, and that this area is effectively 
an extension of the existing channel network. After the 
storm has passed, the channel system eventually reverts to 
its original, perennial extent. Hewlett and Hibbert (1967) 
referred to this idea of a shrinking and expanding source 
area as the variable source area concept. 

Betson (1964) developed a nonlinear mathematical 
infiltration capacity model to analytically equate the 
difference between rainfall and runoff to hydrologic 
variables. Discrepancies between observed and predicted 
results indicated that some important factor had been 
omitted during model development. Further investigation 
showed that the problem centred around the assumption that 
the total watershed area: contributed runoff. The observed 
volume of runoff from the watershed had been converted to 
watershed inches using total drainage area. When adjustments 
were made for smaller contributing, or partial, areas, the 
rainfall-runoff relation occurred according to the 
infiltration capacity concept. Subsequently, it was 
determined that contributing areas for different watersheds 
Supporting various cover types ranged from 4.6 to 86 
percent. The "partial areas" idea helps to explain why 


infiltration capacity of much of the watershed is seldom 


exceeded during storms (Betson, 1964). 
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Ragan (1968) investigated the concept of partial area 
contribution to storm hydrographs by collecting detailed 
information in the vicinity of a 189 m length of a second 
order stream on a forested watershed. Measured inputs to the 
channel included flow entering the channel over an upstream 
Structure, channel precipitation, baseflow prior to storms, 
flow from seeps, and lateral inflows (overland flow, 
interflow and groundwater flow). Outflow from the channel 
section was measured by a downstream control structure. 

Ragan, too, found that only a small portion of the 
watershed contributed flow to the storm hydrograph. The 
contributing area, which is a function of storm duration and 
intensity was very localized - occupying only 1.2 to 3 
percent of the total watershed area. During periods of low 
rainfall intensity most of the flow came from channel 
precipitation and rain falling on wet areas Surrounding a 
series of seeps. During a period of high rainfall intensity, 
flow occurred as saturated porous media flow within the 
lower zone of the forest litter on hillsides, creating a 
larger contributing area. Interflow through mineral soil did 
Not eoccur. 

The idea that only a small portion of a watershed area 
contributes direct runoff has been supported by results from 
other studies. Dickinson and Whiteley (1970) computed the 
minimum contributing area for several storms on the Blue 
Springs Basin in southern Ontario and obtained values 


ranging from 1 to 50 percent, with the majority being less 
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than 10 percent. 

The two concepts of contributing areas referred to in 
the preceding discussion differ in two ways. First, partial 
areas are considered to be more or less fixed in location on 
the watershed, whereas variable source areas expand and 
contract. Second, partial areas supply water to streams by 
means of overland flow, whereas the extended channels in 
variable source areas are believed to be fed by subsurface 
Stormflow (Freeze, 1974). 

Subsurface stormflow which commonly occurs in forested 
areas differs from true groundwater flow in that it flows to 
the stream channel before reaching the general groundwater 
zone. It is more likely to occur where: a) the land is 
Sloping, b) the surface soil is permeable, c) a 
water-impeding layer is near the surface, and d) the soil is 
Saturated (Whipkey, 1965). 

Forest hydrologists (Hewlett, 1961; Hewlett and 
Hibbert, 1963, 1967; Hewlett and Nutter, 1970; Whipkey, 
1965) have studied the subsurface flow mechanism in some 
detail. They believe that the bulk of the average upland 
hydrograph is accounted for through this process (Hewlett 
and Nutter, 1970). To counter the argument that flow 
velocities are too small for subsurface flow to contribute a 
significant volume to direct runoff, Hewlett and Hibbert 
(1967) suggested the concept of translatory flow - a flow 
produced by water displacement. This concept implies that 


"new" water does not have to traverse the soil mass between 
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the point of entry and the stream channel to produce an 
effect on streamflow. Instead, the water induces a "pulse" 
such as occurs when water is added to a soil column drained 
to field capacity. Water flows from the bottom (stream 
channel) almost immediately but it is not the same unit of 
water added at the top (storage). The evolution of the 
variable source area concept and subsurface stormflow theory 
can be traced through a series of papers listed in Hewlett 
(1974) and Whipkey (1965). 

Dunne and Black (1970a,b) examined runoff-producing 
mechanisms on three heavily instrumented hillside plots on 
the Sleepers River watershed, northeastern Vermont. They 
concluded that the major portion of storm runoff is produced 
as overland flow on small saturated areas close to streams - 
a conclusion which supports the partial area concept of 
storm runoff production. Subsurface stormflow was not an 
important contributor to total storm runoff. Freeze (1972b) 
used a deterministic mathematical model to investigate 
mechanisms for generating streamflow. His conclusions, based 
on a number of simulations, were similar to those of Dunne 
and Black (1970a,b). Freeze also specified the conditions 
under which subsurface stormflow becomes large enough to 
contribute significant amounts to storm runoff, namely: 
convex slopes feeding steeply incised channels, and high 
Soil permeabilities. 

It appears then that there are two distinct schools of 


thought regarding streamflow generation in upstream source 
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areas of vegetated watersheds: one supporting the theory of 
overland flow from near-channel partial areas, and the other 
advocating the variable source area concept coupled with 
subsurface flow. However, as Hewlett (1974) points out, if 
one makes allowances for differences in concepts and 
terminology, the two ideas are essentially one and the same. 
There is ample evidence in the more recent literature 
that some or all of the mechanisms described previously may 
be important in runoff generation. The relative dominance of 
_each process is governed by both meteorological and 
watershed characteristics. Hortonian overland flow, for 
example, was the principal mechanism operating on a tropical 
rainforest watershed in Queensland, Australia (Bonell and 
Gilmour, 1978). High intensity rainfalls frequently exceeded 
the average resultant saturated hydraulic conductivities 
below 20 cm depth, which resulted in a perched water table 
and subsurface flow (including "pipe flow") within the top 
20 cm of soil. Widespread saturation overland flow occurred 
when additional rain caused the perched water table to 
emerge at the soil surface. In this case the variable source 
area concept is not applicable (Bonell and Gilmour, 1978). 
Certain soil conditions that prevail prior to storms 
and snowmelt may also be conducive to generation of 
extensive overland flow. During a field study in Yorkshire, 
England, Beven (1978) found that surface runoff occurred on 


large areas of poorly drained soils (which soon reached 


Saturation) during storms. 
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In Canada, snow and ground frost are usually important 
factors governing runoff. Dunne et al. (1976) determined the 
runoff from hillside plots under boreal forest and on tundra 
in the Labrador subarctic, near Schefferville, Quebec. They 
found no evidence of subsurface flow, in the trenches 
prepared for this purpose, during the snowmelt period. The 
soil had been rendered impermeable by concrete frost, and 
runoff took place through the snowpack - primarily through a 
Saturated layer of snow at the base of the pack. 

In contrast to studies which show overland flow to be 
the primary source of storm discharge, there are other 
Studies which indicate that subsurface flow is the only 
mechanism producing runoff (Beasley, 1976; Harr, 1977; 
Megahan, 1972; and Mosley, 1979). Beasley (1976) monitored 
115 storms in northern Mississippi over a three year period 
and found that they produced little overland flow and 
Subsurface flow above the B horizon. Subsurface flow over an. 
impermeable clay layer was as much as 90 percent of rainfall 
during a calendar quarter. It usually began shortly after 
rainfall commenced, even when there was neither saturation 
at the point of outflow nor high antecedent soil moisture. 
Thus the translatory flow concept of Hewlett and Hibbert 
(1967) could not explain this quick response. Instead, the 
response is attributed to the presence of interconnected 
channels through the soil formed by decayed roots and animal 
burrows. This idea is shared by several other workers 


investigating subsurface flow phenomena (De Vries and Chow, 
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1978; Mosley, 1979; and Whipkey, 1965). 

Mosley (1979) used dye tracer experiments on a forested 
watershed in New Zealand to show that water may move through 
Macropores (root channels, pipes) at rates two orders of 
magnitude greater than the saturated hydraulic conductivity 
of the soil. Subsurface flow through the macropores and 
seepage zones in the soil was found to be the predominant 
mechanism of channel stormflow generation. Mosley also found 
that subsurface flow from lower slopes contributed to 
delayed flow. Runoff from partial and variable source areas 
did not contribute significantly to stream discharge, and 
translatory flow apparently was not an important runoff 
generating process on this watershed. 

No overland flow was observed during a study ona 
steep, forested watershed in Oregon (Harr, 1977). Subsurface 
flow and channel interception accounted for 97 percent and 3 
percent respectively of total stormflow. The rapid response 
to precipitation was attributed to translatory flow (Hewlett 
and Hibbert, 1967). Results from stormflows analyses showed 
agreement with the variable source area concept of runoff 
production. 

Weyman (1973) measured saturated and unsaturated 
lateral soil water movement on a convex hillside in 
Somerset, England. He found no overland flow and concluded 
that saturated flow through non-capillary pore spaces 
dominated hillslope hydrographs. Weyman suggested that 


Saturated conditions must be generated within the organic 
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horizons or within the non-capillary spaces of the lower 
soil horizons if lateral soil water flow is to provide storm 
runoff to the stream. Baseflow was dominated by unsaturated 
lateral flow. 

It is evident from the resuits of many hillslope 
hydrology studies (Chamberlin, 1972: De Vries and Chow, 
fe79; Mosley, 1979; Plamondon: et als, 1972; Weyman,« 19733 
and Whipkey, 1965) that macropores Or non-capillary spaces 
in the form of root channels and animal burrows may be 
important factors in storm runoff generation. Plamondon et 
al. (1972) used a water balance procedure to determine 
Saturated flow through macropores of the forest floor and 
soil matrix on a forested site in coastal British Columbia. 
They discovered that usually between 50 and 80 percent of 
the precipitation flowed through the macropores during 
Storms. 

In another experiment on the same watershed De Vries 
and Chow (1978) measured the water potential field during 
wetting and drainage phases of simulated rainfall events. 
Measurements were repeated with the forest floor intact, 
Partly disturbed, and totally removed. It was inferred that 
the flashy response of streams to rainfall events was 
related to rapid subsurface stormflow through root channels. 
Where the forest floor was disturbed there was a tendency 
for less water to move through the root channels, and more 
water through the soil matrix. Consequently, the rate of 


Subsurface flow was diminished. This condition probably 
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resulted from closure of root channel openings during 
disturbance to the forest floor. 

If flow through macropores is indeed an important 
process in storm runoff generation, it is probable that 
Darcy's Law for flow through a porous medium will not apply 
in this case. Instead, the flow may be better described by 
turbulent flow equations. Furthermore, when flow takes place 
through the upper layers, including the litter layer, of the 
forest floor, it also is more likely to occur as turbulent 
rather than porous media flow. Often it will not be readily 
observed because it takes place below the visible coarse 
matrix of the litter, but over the more compact portion of 
the organic horizons. 

In hydrograph analysis groundwater is customarily 
associated with base flow, but there is increasing evidence 
that groundwater may also contribute significantly to storm 
runoff as well. Using environmental tritium, Martinec (1975) 
demonstrated that a substantial part of snowmelt water did 
not leave the basin, although a quantitative water balance 
was maintained. The proportion of subsurface flow 
contributing to snowmelt discharge exceeded 50 percent. 
Apparently snowmelt water infiltrated into the ground and 
"old" subsurface water appeared in the discharge. 

Similar results have been obtained by Sklash and his 
colleagues (Sklash et al., 1976; Sklash and Farvolden, 1979) 
using the Oxygen-18 isotope technique and specific 


conductance methods. They found groundwater to be a major 
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component in virtually all the runoff events they examined. 
Field observations and computer simulations (Sklash and 
Farvolden, 19793) suggest that very large and rapid increases 
in hydraulic head in the near-stream groundwater occur soon 
after the onset of rain. Computer simulations show the 
formation of a groundwater ridge adjacent to the stream in 
response to a rain event with the lag time brief and 
inversely related to the near-stream unsaturated zone 
thickness. 

The same paper also provides a theoretical basis to 
explain the significant role of groundwater in storm and 


Snowmelt runoff generation 


"Along the perimeter of transient and perennial 
discharge areas, the water table and its associated 
capillary fringe lie very close to the ground 
surface. Soon after a rain or snow-melt event 
begins, infiltrating water readily converts the 
near-surface tension-saturated capillary fringe into 
‘a pressure-Saturated zone or groundwater ridge 
(Ragan, 1968). This groundwater ridge not only 
provides the early increased impetus for the 
displacement of the groundwater already ina 
discharge position, but it also results in an 
increase in the size of the groundwater discharge 
area which is essential in producing large 
groundwater contributions to the stream. The 
response of the upland groundwater may become 
important at later times in the runoff event but has 
little influence in the early part of the runoff 
event. 

The groundwater may discharge directly into the 
Stream through the stream bed or it may 1ssue from 
the growing near-stream and/or remote seep areas and 
flow as overland flow to the stream (as in the 
variable source area-overland flow theory). 
Following periods of drought during which the water 
table has fallen far below ground surface, intense 
storms may result in surface saturation from above 
and rain-like overland flow (partial area-overland 
flow) before the water table can emerge." 
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Results from these groundwater studies lend credence to 
the translatory flow concept of Hewlett and Hibbert (1967). 

As new developments of the type described above unfold, 
it becomes necessary to make provision for them in field 
studies and hydrological analyses. For example, the variable 
source area concept is now widely accepted as a valid 
mechanism for generating streamflow. However, techniques to 
identify, map, and predict variable source areas in the 
field on a routine basis have yet to be developed. 

This problem is receiving increasing attention as is 
evidenced from papers by Anderson and Burt (1978), Beven 
(1978), and Dunne et al. (1975). The first of these papers 
(Anderson and Burt, 1978) describes a field study in which 
an automatic tensiometer system was used to monitor soil 
water potential on a single hillside hollow - a location 
usually identified as a variable source area - and on an 
adjacent ridge. Soil water potential maps and flownets 
constructed from field data showed that, after a storm, 
convergent flow into the hollow took place. They also 
established that hollows are the major sources of slope 
discharge while ridge zones were much less important 
contributors. 

While investigating comparative contributions of 
Sideslope and headwater areas to stream discharge Beven 
(1978) also found that convergent hollows are important 
sources of streamflow. These areas of surface soil 


Saturation were associated with ephemeral channels expanding 
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as variable contributing areas of surface flow during storm 

rainfall. The apparent saturated areas were mapped by direct 

observation in the field. Results from the study showed that 
the headwater areas will usually provide higher and earlier 
peak flows per unit area, and more total storm discharge 
than the sideslope areas. 

Dunne et al. (1975) consider repeated field mapping to 
be the best method for determining the size, location, and 
variation of the saturated zone both during and between 
storms. This method should certainly be satisfactory for 
rainstorms but will be impractical where extensive snowpacks 
serve as the water supply. Other methods which can be used 
include: 

1. mapping the saturated areas on the basis of soil 
colouring where the colour is indicative of waterlogged 
soils; 

2. relating saturated areas to areas of "poorly drained" 
soils as shown on soil maps; 

3. using plants as indicators of soil wetness; 

4, relating baseflow, antecendent precipitation index, or 
water table elevation to mapped extent of saturated 


arnea. 


The authors, describing areas for future research, stress 
the need for more field observations relating subsurface 
Stormflow, return flow, and direct precipitation on 


Saturated areas. The process of, and storage opportunities 
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for, saturation overland flow in runoff producing areas also 
needs to be investigated. Techniques are also required for 
estimating the seasonal variation of saturated areas. 

To conclude this section on the watershed system, it is 
worth noting that several investigators who featured 
prominently in the development of the concepts outlined 
above have contributed to a book on hillslope hydrology 
(Kirkby, 1978). This book provides a most comprehensive 
description of the "state-of-the-art" of small watershed 


hydrology. 


D. Watershed And Hydrologic Modelling 

Comprehensive field studies of watershed systems can be 
expensive and time-consuming. Consequently, hydrologic and 
watershed models are used to overcome these limitations and 
to supplement field studies. A watershed model is really one 
type or class of model within the spectrum of a multitude of 
hydrologic models. lene (1973) defines a model as a 
simplified representation of a complex system being either 
physical, analog, or mathematical. Along with the 
development of hydrologic models, a new vocabulary has 
evolved. In their reviews of methods used in hydrologic 
investigations, Amorocho and Hart (1964) and Clarke (1973) 
clarified the definitions of many terms used in hydrologic 
modelling and data processing. The definitions which follow 


are taken from these two papers. 
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Hydrologic studies may be grouped under two areas of 
research: physical hydrology and systems investigations. 
Physical hydrology is concerned with physical processes such 
as those involved in water or energy transfer, whereas 
systems investigations deal with establishing relationships 
between variables, without concern for the physical 
mechanism involved. 

Systems investigations fall into two categories called 
parametric hydrology and stochastic hydrology. Parametric 
hydrology is the development of relationships among physical 
parameters involved in hydrologic events and the use of 
these relationships to generate, or synthesize, non-recorded 
hydrologic sequences. Examples of methods used in parametric 
hydrology include correlation analysis, partial system 
synthesis with linear analysis, general system synthesis and 
general non-linear analysis. Some knowledge of the physical 
system is required when the first three methods are 
employed. 

Stochastic hydrology is the use of statistical 
characteristics of hydrologic variables to solve hydrologic 
problems. A stochastic system differs from a probabilistic 
System, in that it is time dependent. There are two 
important methodologies in common use in stochastic 
hydrology: Monte Carlo methods, and the methods of 
transition probability developed under the general theory of 
Markov processes. The first group is valid for data which 


are statistically independent, and the second is applicable 
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where the state of a system at a particular time depends on 
previous states. 

The different methods of systems investigations share 
two things in common, a) their dependence on historical 
records of parameter values, and b) the assumption of time 
invariance (stationarity) of the hydrologic systems. Both of 
these characteristics place definite theoretical limitations 
on the generality of the solutions. 

Parametric hydrology is usually associated with system 
analyses or "black box" methods in which relationships are 
established between measured input and measured output by 
mathematical procedures. No attempt is made to describe the 
internal mechanisms of the system. In system synthesis, on 
the other hand, the investigator attempts to describe the 
operation of the system by a linkage or combination of 
components, whose presence is presumed to exist in the 
System and whose functions are known and predictable. The 
linkage of components must be made in such a manner that the 
correct output is produced whenever a specified input is | 
applied. System synthesis falls within the realm of physical 
hydrology. 

The nature of hydrologic or watershed models can be 
described in several ways. They may be "deterministic" or 
"stochastic", "conceptual" or "empirical", "lumped" or 
"distributed", and "linear" or "non-linear". If a model or 
Process follows a definite law of certainty but not any law 


of probability, then it is described as aeterminrsticre rt 
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the concept of probability is involved, then the process or 
model re considered to be stochastic or probabilistic. A 
conceptual model is one in which the mathematical functions 
consider the physical processes acting on the input 
variables to produce the output variables. If the model is 
based on obServation and experiment rathern than on 
theoretical considerations, then it is empirical. A lumped 
model is space-independent, whereas a distributed model 
takes into account the spatial variability of input 
variables or model parameters. 

The term linearity may have at least two meanings 
(Clarke, 1973). In the systems theory sense, a model is 
linear if the principle of Superposition holds: that is, 
given that y,(t), y2(t) are the outputs corresponding to 
inputs x,(t), x2(t), a model is linear if the output 
@orresponding itoisinput’ «x ¢(ty)i it ex, (tt) chs: yaGtdeet cy alt) o This 
is the sense in which linearity is most widely used in the 
hydrologic literature. Alternatively, the model may be 
linear in the statistical regression sense if it is linear 
in the parameters to be estimated. 

Improved computer technology and the International 
Hydrological Decade (1965-1974) have provided the stimulus 
for a considerable amount of research in watershed and 
hydrologic modelling. Several mathematical or digital 
watershed models were published in the sixties and 
Seventies, the best-known being the Stanford Watershed Model 


(SWM) (Crawford and Linsley, 1966). A flow chart for the SWM 
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is shown in Fig. 5. 

This model, the parameters of which are optimized, can 
be used to simulate watershed response to weather 
modification and urban development; to design 
hydrometeorological networks; and to predict streamflow from 
ungauged watersheds. Streamflow is calculated at several 
locations in the stream channel called flowpoints. The area 
upstream from each flowpoint is divided into segments so 
that there are one or more segments for each recording rain 
gauge. The model continuously calculates streamflow at each 
flowpoint from rainfall in each successive watershed 
segment, and from flows measured or calculated at upstream 
flowpoints. 

Major data inputs for the Stanford model are 
precipitation and potential evapotranspiration. Further 
meteorological data are required if the snowmelt subroutine 
is used. Available moisture and potential evapotranspiration 
data inputs are used together with cumulative frequency 
distributions of infiltration and evapotranspiration 
opportunity to determine detention, infiltration, 
interflows, and evapotranspiration. 

Three subsurface water storages are defined: the upper 
zone, the lower zone, and groundwater storages. The upper 
and lower zone storages control overland flow, infiltration, 
interflow, and inflow to groundwater storage. The upper zone 
simulates the initial watershed response to rainfall and is 


of major importance for smaller storms and for the first few 
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hours of larger storms. The lower zone controls watershed 
response to major storms by controlling long term 
infiltration. Base flow to stream channels is supplied from 
groundwater storage. None of the soil water storages have 
fixed capacities, and evaporation and transpiration may 
occur from any of them. Streamflow hydrographs are the most 
important outputs from the Stanford model which has been 
tested using data from over 300 basins ranging in size from 
0.65 to 800,000 km? (Linsley, 1976). If reliable input data 
and accurate streamflow records are available, agreement 
between observed and simulated flow is good with annual flow 
volumes within +5 percent and peak flows within +10 percent. 
Another well-known hydrologic model is the Streamflow 
Synthesis and Reservoir Regulation (SSARR) model (Rockwood, 
1958, 1968; US Army Corps of Engineers, 1971). SSARR is a 
deterministic hydrologic mathematical river system model 
designed to simulate flows from large basins. It uses basin 
storage routing procedures to convert precipitation excess 
over the basin to represent streamflow. Soil moisture, 
evapotranspiration, and base flow indices, together with 
empirically derived relationships, are used to determine 
soil moisture increases, soil moisture reductions, and 
baseflow contributions, respectively. The model also 
contains snow accumulation and snowmelt components. 
Although SSARR is designed primarily for streamflow 
forecasting and reservoir regulation it can also be used in 


the development of design floods, or in the extension of 
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missing or broken records of streamflow. Examples of 
important applications of SSARR include simulations of flow 
from the Columbia River Basin in Western North America, and 
from the Lower Mekong River in Southeast Asia. 

In contrast to the engineering approach evident in the 
Stanford and SSARR models, Freeze (1971) developed a 
completely physically based, theoretical hydrologic response 
model and used the entire set of equations of flow through 
heterogeneous, anisotropic porous media to obtain solutions 
to the transient, saturated/unsaturated flow problem for a 
groundwater basin. 

In this model, the watershed or sideslope is defined in 
terms of a block-centred nodal grid (Fig. 6), and mesh 
Spacings which can be varied. Inflows such as rainfall and 
outflows such as evaporation are simulated at specified 
boundary flux nodes. At other boundary nodes, boundary 
conditions may be imposed that specify hydraulic head or 
no-flow conditions. The basic information required to 
operate the model is related primarily to the properties of 
water, soil and the geologic formations, namely: the 
compressibilities of soil and water; porosity, specific 
permeability, density of water, and water content. The last 
four variables are functions of pressure head. 

The fundamental equation of flow for the model is a 
nonlinear parabolic partial differential equation. An 
iterative numerical scheme that employs an implicit finite 


difference formulation, and known as the line successive 
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overrelaxation technique, was used to obtain solutions to 
the equation. Output takes the form of plots of pressure 
head, total head, and moisture content fields for any cross 
section at any time step. Base flow is determined from 
flownet analysis. 

Freeze (1972a) subsequently coupled his subsurface flow 
model with a one-dimensional, gradually varied, unsteady 
channel flow model to study the mechanism of base flow 
generation and the nature of watershed response in streams 
dominated by base flow. The full set of nonlinear, 
hyperbolic partial differential equations for shallow water 
flow were solved uSing an explicit finite difference 
technique known as the single step Lax-Wendroff method. The 
combined model operates (Freeze, 1972a) as follows: 

"...internal coupling is carried out by using the 
Stream depth as the convergence quantity. The 
subsurface model is solved at each time step by 
uSing the stream depth from the previous time step 
as the specified head condition at the stream 
boundary. The calculated outflow from the subsurface 
System becomes the lateral inflow to the streamflow 
model for that time step, and a new stream depth 
profile is calculated. This new profile replaces the 
old to set the heads for a solution to second cycle 
subsurface. This alternating cycle is continued 
until the stream depth and the specified head 


boundary values converge to within a predetermined 
tolerance." 


The uncoupled subsurface flow model was designed 
primarily as a research tool and has been applied to many 
hypothetical hydrologic situations. It has also been applied 


to the real world conditions of Reynolds Creek watershed, 
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Idaho (Stephenson and Freeze, 1974). Calibration of the 
model consisted of trial-and-error matching of computer 
output with field measurements of check parameters, during a 
series of simulations in which some of the input parameters 
were varied. Results from the field study indicate that, 
within the limitations of the model and data availability, 
the model adequately simulates natural conditions. However, 
because of: its extreme complexity, its requirements for 
large amounts of data, limitations in theoretical 
development, constraints on calibration procedures, and 
computer storage limitations, it cannot be considered 
Suitable for routine operational use. 

The interaction between biologically controlled and 
physically controlled hydrologic processes, particularly in 
interception, throughfall, evapotranspiration, snow 
accumulation, snowmelt, and soil water redistribution 
Simulations, features prominently in some watershed models. 
Such models are most useful for investigating plant-water 
relationships, and for simulating the effects of land use on 
hydrologic processes for watersheds where vegetation is an 
important consideration. 

The United States Department of Agriculture Hydrograph 
Laboratory model of watershed hydrology, USDAHL~74 (Holtan 


1., 1975) was designed to serve the needs of 


et 
agricultural watershed engineering. It employs plant growth 
indices in both evapotranspiration and infiltration 


simulations, and requires more detailed information on soil 
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physical properties - particularly in the plant root zone - 
than does either Stanford or SSARR models. Soils on each 
watershed are grouped by land capability classes to form 
hydrologic response zones for computing infiltration, 
evapotranspiration, and overland flow. Three zones are 
identified on a watershed: uplands, hillsides, and 
bottomlands, therefore the model can be categorized as a 
distributed one. Computations are based on the assumption 
that some of the runoff will cascade over successive zones. 

The program calculates the extractions from a soil 
layer through evapotranspiration and free drainage, and then 
records the resulting increase in available storage during a 
given time interval. Reductions in available storage are 
computed when infiltration occurs from surface depression 
Storage, or following rainfall or overland flow. Water 
inputs applied at rates in excess of the infiltration rate 
are routed as overland flow. Channel flows and subsurface 
return flows are aamhace by solving the continuity equation 
and a storage function using recession flow analysis. 

A very detailed biophysical hydrologic model has been 
assembled by Huff and his colleagues from the Oak Ridge 
National Laboratory, Tennessee (Huff et al., 1977). It is 
called the terrestrial ecosystem hydrology model (TEHM), and 
is essentially a combination of two models that were 
developed independently: the Wisconsin Hydrologic Transport 
Model and an atmosphere-soil-plant water flow model 


(PROSPER). The first is based in part on the Stanford model. 
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In PROSPER, a water balance is applied to a stand of 
vegetation on a layered soil. The model, which is 
one-dimensional, uses a combined energy balance-aerodynamic 
method to derive an equation for evapotranspiration as a 
function of a resistance to vapor transfer which is 
characteristic of evapotranspiration surfaces. Because of 
the close parallel between the equations of flow for 
electricity and those of water flux through a 
soil-plant-atmosphere continuum, an electrical circuit 
analogy (Fig. 7) was developed so that standard methods for 
solving electrical circuit problems could be applied. 
PROSPER has been used as a completely self-contained 
model to simulate evapotranspiration and annual drainage 
from a mature hardwood forest at Coweeta in the southern 
Appalachians (Swift et al., 1975). Because the movement of 
water between 90 cm depth and the stream is not simulated in 
PROSPER, the drainage term was equated to streamflow for 
annual periods. For the hardwood forest, simulated annual 
drainage was within 1.5 percent of measured streamflow. 
Simulations were also performed for a 16-year old white pine 
plantation and a 1-year old clearcut area. The model 
estimated a reduction in drainage of 20 cm for the 
plantation and an increase of 36 cm for the clearcut. These 
values correspond to measured annual streamflow changes of 
-20 and +38 cm respectively. The results from these 
simulations also indicated that summer evapotranspiration 


for hardwood and pine were identical, and that during winter 
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Figure 7. Schematic of PROSPER (After Huff et al., 2972). 
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and early spring, losses due to evapotranspiration were 
greater for pine. Simulations showed that evapotranspiration 
was considerably less, and soils much wetter, on the 
clearcut areas than on the forested areas. 

The subsurface component of the TEHM takes into 
consideration the variable source area concept, by relating 
subsurface drainage rate and variable source area runoff. As 
Huff et al. (1977) explain, when the root (or surface) zone 
drainage rate exceeds a certain threshold value, the 
fraction of source areas draining to the channel increases 
linearly with drainage rate until a specified upper limit is 
reached. This upper limit defines the full extent of 
variable eontce areas. The lower limit represents the 
permanent part of the source areas that are linked to the 
Stream in the basin, and will contribute subsurface drainage 
through the layers below the root zone at all times. Beneath 
the zones simulated in PROSPER, soil water transmission 
zones are defined in which soil hydraulic conductivity is 


related to soil water content by the expression 


K(6)=KS , ei? — 4) 1<is3 (2.1) 
where 
K(6) hydraulic conductivity at water content 6, 
KS | Saturated conductivity parameter, 
aj curve fitting parameter, 
6 water content, 
Oo, water content at saturation, 
1 section number of portion of K vs. @ curve. 
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Hydrograph separation techniques are used to determine 
outflows from groundwater storage. 

The channel flow portion of the model combines 
Simulated flow components such as direct channel 
precipitation, impervious area runoff, runoff from permanent 
and variable source areas, and groundwater inflow from each 
portion of the watershed. It then generates hydrographs at 
specified points in the channel system. 

Output from the TEHM consists of monthly and annual 
water balance summaries, and daily values of selected 
hydrologic variables for each watershed segment. Water flux 
terms such as infiltration, direct runoff, drainage to 
deeper soils, evapotranspiration, lateral subsurface flow, 
and net flux from the second to third soil layer, as 
calculated by PROSPER for each day are tabulated. 

The monthly summary includes a "deep soil balance” term 
in addition to the terms which normally appear in water 
budget equations. This term is related to storage changes in 
the soil water transmission layers. 

Four summaries appear in the annual summary: runoff, 
precipitation/evapotranspiration, moisture status, and water 
budget. The moisture status summary is intended as an index 
to total water content of a vertical column extending from 
Canopy to bedrock, but it is more likely indicative of 
overall water storage in the unsaturated soil column (Huff 
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The TEHM has been tested using data from the Walker 
Branch ne renateam a 97.5 ha drainage area located in the 
Ridge and Valley province of east Tennessee. A predominantly 
oak-hickory forest covers the basin. During simulations the 
model generated hydrographs for stormflow, base flow and 
total flow. The total flow data were compared with the 
observed hydrographs, while stormflow and baseflow were 
compared with quickflow and delayed flow, respectively. 
Quickflow and delayed flow were calculated by the method of 
Hewlett and Hibbert (1967). In all cases good agreement was 
obtained between simulated and observed results. Monthly and 
annual summaries were also provided by the model. 

A special feature of the terrestrial ecosystem 
hydrology model is that it can be used in environmental 
chemistry studies. For this purpose TEHM is combined with 
Subroutines for atmospheric transport, soil chemistry, 
exchange of heavy metals, a plant-growth model, and a 
mineral uptake model. 

A number of hydrologic or watershed models have been 
developed specifically to meet the needs of forest 
hydrologists or watershed managers. Multiple regression and 
covariance analyses are perhaps the earliest and simplest 
examples of such models. Anderson (1960) used regression 
models to relate forest shade and radiation to snow 
accumulation and melt. The results indicated that shade 
could be maximized and back radiation minimized by cutting 


east-west strips in the forest, with successive cuts 
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proceeding toward the south. This idealized forest stand 
structure, referred to by Anderson as a wall-and-step forest 
(Fig. 8), would maximize both snow accumulation and delay in 
snowmelt, and therefore would be favourable for water 
production. 

Leaf and Brink (1973a) developed a snowmelt simulation 
model for Colorado subalpine watersheds and used it (Leaf 
and Brink, 1972, 1973a) to simulate the probable Seeecre of 
timber harvesting on snowmelt in the subalpine forest. The 
model was subsequently expanded (Leaf and Brink, 1973b) to 
Simulate the total water balance on a continuous, year-round 
basis. It was designed to simulate watershed management 
practices and their effects on hydrologic systems. 

For simulation purposes, the study basin is divided 
into several hydrologic response units. Hydrologic responses 
are computed for each unit, then weighted according to their 
respective areas and combined to produce a composite record 
of hydrologic system behaviour on a watershed basis. Input 
to the response system is derived from snowmelt and 
rainfall. Once evapotranspiration requirements have been 
Satisfied, additional input satisfies soil mantle recharge 
requirements. When field capacity is reached, residual input 
becomes water available for streamflows. Excess water is not 
routed through the soil mantle. 

For this model, the structure and properties of the 
forest stand play an important role in the simulation 


process (Leaf and Brink, 1973b). Forest cover density, for 
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PROGRESSIVE CUTS 
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Figure 8. Wall-and-step forest - a cutting pattern for su ele deat 


yield from snow zone forests (After Anderson, 1960). 


We 


tee tas 
sviaesenioed 


2 frau hs 
(Oat) | 
iti bik 
fete - 
ee ae a. 
ie ik es 
Y ad BAS ‘ 
Bi Fre 
he rit 
eb 
gw : 


74 


example, is used to determine the amount of snow intercepted 
by the forest canopy and vaporized from it. It is also 
required to compute evaporation from the snowpack beneath 
the trees. Reflectivity of the forest stand is a function of 
forest density, and is used to calculate evapotranspiration 
from forest and open areas. Cut blocks and other clearings 
are simulated by setting the forest cover density equal to 
zero. 

The Colorado subalpine hydrologic simulation or water 
balance model serves as the core system of the Land Use 
Model also developed by Leaf and Brink (1975). The latter 
model was designed to extend the capabilities of the former 
Over a much longer time period - from a few years to 120 
years or more (the rotation age of subalpine forests). The 
package provides greater flexibility for simulating a 
variety of timber harvesting and weather modification 
combinations. 

The interactions between timber management and 
watershed management in the subalpine forests of Colorado 
have also been simulated (Leaf and Alexander, 1975). This 
was accomplished by linking timber models for lodgepole pine 
and spruce-fir forests to the water balance and land use 
models. Linkage was achieved through the forest cover 
density variable. Projected changes in water yield were 
obtained for several different forest management practices. 
The investigators concluded that this method for projecting 


long-term yields of both wood and water shows that the 
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Strategies for optimum water production are compatible with 
the conversion of old growth forests to managed stands for 
timber production. The major disadvantage of both the land 
use and the timber/water model is the inability to verify 
the long-term projection during the short term. 

Another model designed to predict long-term timber and 
water yields from snow pack zones of lodgepole pine 
subalpine watersheds along the Colorado Rocky Mountain Front 
Range was developed by Betters (1975). The model operates on 
an annual time scale and incorporates several empirical 
relationships to simulate forest growth and water production 
on a watershed. Precipitation, potential and actual forest 
evapotranspiration, soil water holding capacity, and 
proposed timber-harvesting schedule are the most important 
variables and parameters used by the hydrologic submodel. 

Past records of timber inventory and runoff values for 
Little Beaver Creek and Fool Creek were used to test the 
model. Both watersheds have long-term streamflow records. 
Results from the tests indicated that the predicted timber 
inventory estimate was 8 percent higher than the actual 
timber inventory estimate. On Beaver Creek, the observed 
water yield exceeded the predicted yield by 13 percent. Two 
time periods were simulated for Fool Creek ~- before harvest 
and after harvest. Predicted water yields exceeded observed 
by 16 and 10 percent, respectively. 

Simulations showed that on a per hectare basis, after 


harvest the water yield from the uncut portion of a subunit 
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decreases. The net result, however, is an increased yield 
because the large increase in snowpack combined with reduced 
evapotranspiration on the cut area exceeds the loss 
Sustained on the uncut portion (Table 5). 

To date there have been few attempts to simulate the 
effects of forest succession on water yield. Forest 
succession is the gradual replacement of one forest tree 
Species or plant community by another. Jaynes (1978) 
developed a lumped, deterministic model (ASPCON) to simulate 
the hydrologic impacts of forest succession in the Rocky 
Mountain area of Utah. The model calculates weekly water 
budgets throughout one water-year, using inputs of 
precipitation and air temperature. The model predicted 
(Table 6) a net loss of 8.6 cm in streamflow when aspen 
assumed dominance over grass-forb’* cover. An additional loss 
in streamflow of 11.7 cm was predicted when conifers 
replaced the aspen. 

Attempts to incorporate relatively new hydrological 
concepts into hydrologic or watershed models are exemplified 
by the work of Beven and Kirkby (1979). They formulated a 
physically based, deterministic model that provided for 
Simulation of variable source areas. The model was designed 
to predict the hydrologic response of ungauged basins 
without resorting to optimization procedures. Only 
information derived directly from the basin itself was used. 


Application of the model to an 8 km’ watershed showed that 


*An herb other than grass 
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there was fairly close agreement between simulation output 
and observed values. 

The variable source area concept has been incorporated 
into other hydrologic models - models designed specifically 
to evaluate the effects of forest land use on the hydrology 
of an area (Troendle, 1979; Chanasyk, 1980). 

Troendle's (1979) simulator contains three main flow 
components: subsurface flow, channel interception, and 
impervious area overland flow. The simulation watershed is 
partitioned into a number of subbasins and each subbasin is 
subdivided into segments and horizons. The segments run from 
ridge-top to stream channel, and the horizons correspond to 
soil layers. After the system has been initialized, an 
explicit finite difference scheme is used to determine the 
volumetric water content of each element, at regular time 
intervals. During each iteration, a test is made for 
Saturation in the surface layer beginning with the first 
Segment next to the stream channel, and then proceeding 
upslope until an unsaturated element is detected. At this 
point, the saturated surface is considered to be part of the 
expanding channel system, and the remaining unsaturated 
portion of the slope is partitioned again for the next 
iteration. A reverse procedure is used during storm 
recession simulation, when channel contraction takes place. 
Tests are made on previously saturated zones beginning with 
tridge-top elements and proceeding towards the channel. If an 


element is unsaturated, it is added to the unsaturated slope 
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which is re-partitioned for the next iteration. This concept 
of expanding and contracting channels is the central idea 
behind the variable source area model. 

Troendle (1979) compared results from his simulator 
with data from a 38 ha watershed, in north central West 
Virginia, without resorting to calibration. He found that 
Simulated storm response varied from 45 to 125 percent, and 
averaged 89 percent of one streamflow. The simulated time to 
peak was about 4 percent below the observed time to peak. 
More than 90 percent of the simulated streamflow was of 
Subsurface origin. 

The idea of expanding and contracting stream channel 
Systems for low order watersheds is also incorporated into 
Chanasyk's (1980) land use hydrologic model. The model, 
called SLUICES (Soils and Land Use affecting Interflow and 
Creating Effects on Streamflow), uses a square element grid 
System, together with seven easily obtainable or optimized 
parameters. The most sensitive of these parameters are: 
hydraulic conductivity, saturation moisture content, field 
Capacity, and drainage coefficient. The model is also 
ate to antecedent soil moisture conditions and to 
evapotranspiration. Model operation is based on the 
existence of three storages: unsaturated flow storage, 
Saturated flow storage, and overland flow storage. When the 
unsaturated storage is completely filled, it becomes 
Saturated storage. Additional water creates overland flow 
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Chanasyk (1980) applied his model to the Jamieson Creek 
watershed, located near Vancouver, British Columbia. 
Agreement between recorded and simulated discharge 
hydrographs for seven storms ranged from very good to fair - 
depending on whether the storms were preceded by wet or dry 
antecedent conditions , respectively. He found that 
streamflow can be augmented by logging , and that the 
increase waS proportional to the area clearcut. 

It is evident from the preceding discussion that 
hydrologic and watershed modelling has advanced considerably 
during the past two decades. Extremely sophisticated 
mathematical models are currently available that can be used 
to simulate forest and watershed management situations. It 
is also clear that the simulation of certain processes, 
notably the subsurface flow components, needs to be improved 
if the goal of creating a completely physically based 
hydrologic or watershed model is to be achieved. Groundwater 
hydrologists, petroleum engineers, and soil physicists have 
paved the way to meeting this objective by deriving and 
solving the complicated mathematical expressions which 


describe subsurface flow. 
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E. Development Of Subsurface Flow Equations And Methods For 


Obtaining Their Solutions 


Saturated flow 

The equation for groundwater flow developed by Darcy 
(1856) from his experiments on water flow through filter 
sands may be considered as the forerunner of most of the 
Subsurface flow equations developed subsequently. It can be 


written (Davis and DeWiest, 1966) as 


Q=KA(H,-H2 ) =-KAdGH G22) 
al dl 

where 

Q flow rate, 

K hydraulic conductivity, 

A cross sectional area of flow, ; 

H; energy per unit weight of fluid or hydraulic head; 
in the case of water = Z+Ptarbitrary constant, 

Y 

Z elevation above an arbitrary datum plane, 

P pressure sustained by the fluid in the pores of 
the medium, 

7 Specific weight of the fluid, 

agH/dl hydraulic gradient. 


The Darcy equation is used together with the equation 
of continuity (Jacob, 1949) to obtain Laplace's equation for 


Steady, incompressible flow 


a2H+02H+92H=0 (2.3) 
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where H is the hydraulic head. 
This is an elliptic, second-order partial differential 
equation. 

Fifty years after Darcy's Law was first formulated 
Thiem (1906) used Darcy's Law in terms of polar coordinates 
to determine the cone of depression for groundwater in the 
vicinity of a discharging well. However, it could be used 
only when steady shape or equilibrium conditions prevailed. 

Later, Theis (1935) used the analogy between heat flow 
and groundwater flow to develop an equation applicable under 
nonsteady or nonequilibrium conditions. The time variable is 
introduced into the equation, and an analogy is drawn 
between sources and recharging wells, and between sinks and 
discharging wells. The Theis equation can be used to 
determine drawdown before equilibrium conditions occur. The 
Same fundamental differential equation and its solution were 
obtained by Jacob (1940) from Darcy's Law and the continuity 
equation. Jacob (1949) subsequently derived the 
three-dimensional form of the general differential equation 


for unsteady groundwater flow 


07H+07H+07H=S0H=pq(nf+a) dH (2.4) 
Kot K 
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where 

H hydraulic head, 

S Storage coefficient, 

K hye@raulic conductivity, 

p Fluid density, 

g gravitational constant, 

n porosity, 

B GCompressibility of fluid, 
a compressibility of medium. 


This is a parabolic partial differential equation. Cooper 
(1966) was able to verify Jacob's derivation and resolve 
certain contradictions by considering the Z (vertical) 
coordinate as a deforming coordinate. 

Before the introduction of sophisticated computers, 
equations (2.3) and (2.4) or their polar coordinate 
equivalents were solved analytically by the methods 
developed by Thiem (1906) and Theis (1935). However, since 
the mid-fifties groundwater hydrologists have been turning 
increasingly to numerical methods for solutions to these 
equations. Freeze and Witherspoon (1966) used both 
mathematical and finite difference techniques to solve the 
Laplace equation [equation (2.3)] for a three-dimensional, 
nonhomogeneous, anisotropic regional groundwater basin. They 
compared the method of analytical separation of variables 
and Fourier Series technique employed in partial 
differential equations theory with a numerical finite 
difference approach, and found the numerical method more 
versatile, mathematically simpler, and better suited to 


computer-oriented methods of data storage. 
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Laplace's equation is the governing equation in the 
problem of steady seepage through a homogeneous earth dam 
and as such has been solved by means of finite difference 
and relaxation methods (Finnemore and Perry, 1968). 
Laplace's equation and a finite difference technique have 
also been used to predict the effect of a proposed reservoir 
on aquifer water levels (Remson et al., 1965). 

Equation (2.4) also has been solved using numerical 
methods. Pinder and Bredehoeft (1968) used an 
alternating-direction implicit technique to solve the 
transient flow problem in the case of an aquifer in Nova 
Scotia. Their calculated values compared favorably with the 
analytical solutions for homogeneous, isotropic aquifers of 
Simple geometry. The same iterative technique has been used 
to solve other practical (Trescott et al., 1970) and 


theoretical (Bredehoeft and Pinder, 1970) problems. 


Unsaturated flow 

It will be observed that no attempt has been made to 
define the saturated and unsaturated zones or the water 
table. The omission is deliberate inasmuch as we wish to 
pursue a continuum approach to the subsurface flow problem. 
However, it is essential to differentiate between saturated 
and unsaturated flow, and to discuss factors related to this 
difference. This can be achieved best by first discussing 


the concept of potential. 
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For the equations discussed so far(2.2 to 2.4), the 
term 'H' is referred to as hydraulic head. It is also called 
Ppiezometric head. Because this term has the dimension of 
energy per unit weight, H may be regarded as the potential 
energy of water at a point (Jacob, 1949). According to 


Remson et al. (1971) 


",..The [total] potential at a point is a measure of 
the work required to transport a test body of water 
to that point from a specific reference state ina 
soil-water system which is in a state of rest. When 
described in terms of energy per unit weight of 
water, potential has the dimension length and is 
referred to as "head". The gradient of total 
potential is proportional to the water-moving 

BOECES sksce. 


Evidently, the terms: hydraulic head, piezometric head, 
total potential, and piezometric potential can be used 
interchangeably. This concept is applicable to both 
Saturated and unsaturated flow. 

The total potential (®@) at a point can be written 


(Remson et al, 1971) as 


B=Y, +V, HV, tVe tH, tHe WA) 
where 
Wg gravitational potential, 
Yp hydrostatic pressure potential, 
Wo osmotic potential, 
Wo adsorption potential, 
VY; thermal potential, 


We chemical potential. 
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The gravitational potential is related to the position of 
the point and depends only upon its height above datum. Thus 
Mepis 1dentical’ to the’ elevation head (z) of that point. 

The hydrostatic pressure potential is considered to be 
zero at the water table where pressure is atmospheric. 
Therefore all pressures in subsurface water are gauge 
pressures. Below the water table at static equilibrium, 
hydrostatic pressure potential increases with increasing 
depth. In the case of unsaturated media, the hydrostatic 
pressure potential is negative (i.e. it represents a 
Suction) because a suction must be applied in order to 
withdraw water °. 

The osmotic potential is related to the effect of two 
solutions of unequal concentrations being in contact through 
a semipermeable membrane. Water moves through the membrane 
from the less concentrated to the more concentrated 
solution, thereby increasing the pressure in the more 
concentrated solution. This pressure is termed the osmotic 
pressure potential. 

Frequently, cations become dissociated from soil 
particles of clay size leaving the surface of the soil 
particle with a negative charge. The attraction of particles 
having opposite electrical charge results in an increased 
concentration of cations in the vicinity of the negatively 


charged particles. It is this increase in cations which 


‘This concept of treating suction as a negative hydrostatic 
potential for unsaturated conditions will be maintained 1n 
this study. 
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results in an osmotic force moving water into the space 
between the particles. The osmotic pressure potential is 
important in soils that shrink upon drying. 

The adsorption potential is due to the attractive 
(adsorptive) forces between the matrix water and the matrix 
particles and are prominent only when the soil becomes 
extremely dry. 

Thermal potential has appreciable effect only in the 
case of vapour diffusion or capillarity. 

Chemical potential is due to the osmotic energy of ions 
free in the aqueous solution. This is distinct from the 
potential due to the osmotic energy described above. 

If we assume isothermal conditions and uniform solute 


concentration, the total potential at a point becomes 
B=Y,+Vy ty ta (2e6) 


In the case of saturated flow the last two terms on the 
right hand side (W,. and Wa) are disregarded. For unsaturated 
flows, all the terms are retained. However, for 
soil-moisture ranges usually encountered in the field, Wa 
can be neglected. In soils that shrink on drying Wo assumes 
greater significance than Wp; in nonshrinking soils the 
reverse is true. Generally speaking, the last three terms on 
the right hand side of equation (2.6) cannot be evaluated 
Separately, so commonly they are lumped together and 


referred to collectively as the "capillary" potential. The 
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basic idea of capillary potential which, together with 

"capillary conductivity," forms the basis of modern soil 

water physics was first set forth by Buckingham (1907). 
For the remainder of this treatise, the following 


definition of total potential (#@), and notation, will be 


used 

G=y+z (2879) 
where 
W hydrostatic pressure potential (or pressure head), 
z gravitational potential (or elevation head). 


The extension of Darcy's Law to unsaturated flow was 
first proposed by Richards (1931) who assumed that the 
conductivity K, as well as the water content could be 
treated as (non—hysteretic) functions of pressure potential 
(Wy). If the soil water is assumed incompressible, the 
general equation for unsaturated, unsteady flow can be 


Written as 


0 (K(w)ab]+a_ [K(w) ob]+0_[K(y) ab]=90 CZ otshy 
Ox ox dy Ch or Z t 
where 
K(y) hydraulic conductivity as a function of 
pressure potential (Ww), 
$ total potential, 


6 volumetric soil water content. 
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This equation, sometimes referred to as the Richards 
equation, is highly nonlinear. Since most of the solutions 
for unsaturated flow problems have been obtained using the 
one-dimensional case, it is useful to write the equation for 


vertical soil water movement 


8 _[K(y) a@]=00 (2.9) 
dz or rot 

Since = W+z (equation 2.7), we can write 
d_[K(y)3 (wtz) ]=200 (2.10) 
dz OZ ot 

or 
a_[K(W) dwt+K (wy) ]=96 | (254) 
dz dz dt 


Childs (1936) hypothesized that water movement in the 
unsaturated zone is determined by the moisture concentration 
gradients, which implies that water moves according to a 


diffusion equation 


dc=-kd7c (C2Rai2) 
at aax- 

where 

C moisture concentration, 

k diffusion coefficient (constant), 


0 weight of dry matter per unit volume. 
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When ¥ and K are single-valued functions of 6, equation 


(2.11) becomes 


d_[Dd6+K(6) ]=96 (27513) 
OZ Z c 


where D=K(@)dW/00 or moisture diffusivity (Childs and 
Collis-George, 1950). 
The term diffusivity and the symbol D are used because the 
form of the equation is the same as that of Fick's Law of 
diffusion; there is no implication that molecular diffusion 
1s or iS not involved as a mechanism (Miller and Klute, 
1967). Equation (2.13) is commonly used as an alternative to 
equation (2.11). 

In order to introduce one further concept related to 
unsaturated flow, the right hand side of the foregoing 


equations will be written as 


(214) 


Qa 


262 
ow 


— 


Q 


t 


The term 06/dw is defined as the specific moisture capacity 
C(@) and represents the slope of the water retention or 
moisture characteristic curve (6 vs Ww). 

It is evident from the preceding discussion that the 
main difference between saturated and unsaturated flow is 
the dependence of K, Ww, and C on soil water content in the 


unsaturated case. At saturation, W is no longer a function 
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of 6, K iS a constant, and C=0. 

A basic feature of unsaturated flow problems is that K 
and w are not single-valued functions of @ (although the 
Opposite assumption is often made). This phenomenon is 
termed ‘hysteresis’, and implies that a wetting soil will 
have a characteristic curve different from that of the same 
s0l1l when it is drying. Ree tnernene, there may be any number 
of intermediate or scanning curves which are governed by the 
Soll water content at the time when change (from wetting to 
drying, or from drying to wetting) occurs. Usually the term 
"characteristic curve' is reserved for the limiting wetting 
and drying curves. 

The first Serious attempt to solve the unsaturated flow 
equation appears to be that of Klute (1952). He used an 
iterative procedure to obtain moisture content distributions 
in semi-infinite columns of sand and clay after varying time 
moeervals, in the case of horizontal flow..Philip .(1957) 
also used an iterative method to solve the one-dimensional, 
horizontal flow equation, and described a numerical 
procedure to solve the vertical, unsaturated flow problem. 

Hanks and Bowers (1962) were perhaps the first to 
investigate layered soils using the unsaturated flow 
equation, which they solved by means of the Crank-Nicolson 
implicit finite difference technique. The effects of 
continuous rainfall on soil water content were examined by 
Rubin and Steinhardt (1963). They used a linear 


extrapolation-based linearization of the finite difference 
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equation to solve the Richards equation, and concluded that 
ponding at the surface results only if rain intensity 
exceeds the soil's saturated hydraulic conductivity. Whisler 
and Klute (1965) and Rubin (1967) took hysteresis into 
consideration when evaluating unsaturated flow problems. 
This represented a significant advance in unsaturated flow 
analysis. An explicit-implicit difference scheme was used by 
Wang and Lakshminarayana (1968) to investigate unsteady soil 
water movement in a nonhomogeneous unsaturated soil. 

Freeze (1969) reviewed the more important papers 
dealing with application of numerical methods to the 
unsaturated flow equation, and tabulated the significant 
features of each paper. He presented a mathematical model 
for one-dimensional, vertical, unsteady infiltration or 
evaporation above a recharging or discharging groundwater 
flow sytem. Freeze also used hypothetical cases to examine 
the effects of soil type, rainfall intensity and duration, 
antecedent soil water conditions, groundwater recharge or 
discharge rate, depth to water table, and depth of ponding, 
on the saturated-unsaturated flow system. 

Hanks et al. (1969) used a modification of the 
procedure developed by Hanks and Bowers (1962) to estimate 
one-dimensional infiltration, redistribution, evaporation, 
and drainage of water from soil - taking hysteresis into 
account. They found good agreement between measured and 


Calculated values. 


: r re ah ; on pr 
fet bebulones bas <npisgupe Seaeisin Sas oer 


P) 


_ 
we wt 


~ 


v+iengier agen’ eine ea iirey: eoetada 


‘zenedenn A003 Cee “hitwh Bre cate 


'Jemagran 


Ds ' 2 
1. Tihs teuosas” LG pet LH De ie oe 


Ve noted eeu ees MCitaesti tins 


lagisdedit Bas tied cise mela? SSE 


23283) Satrervodva-HSeu. ict sine ba 


4656 bre elas ladsew, 48 sata abst 


ae 


=e 


+suonoS ottustbed botnet end 


(O$3, Boteauseeng” pny? eu Leres pee, oe 


sf so ts5¢e7e2 B, iP 3: 41gKe8 
stieottesvil Of (e°4)). SASPe Tense 
/ | : 
Ssjayvt seal euosdaouasinds ane 
f5sIOS ion «rd Bevesues (eat 
sbonts seciiemyun bo solesoitegse 


a 


Liens VoRSssenu ,la012Say rants . 7 
ae 


s vtbedeini- Lisintas joe tise » & 


‘ 
TS 7 othe sD pre Qis bia SPIEM 


ri 


Ce ee 
> a 


: a we 
athe, 
a wes): 
smpteve Seebich vei he pti é 338 " 
ry Mee Pie Wel ai ~ 
sitbhom @ fsa (2ae7) eS ra F 


(REEL) 23 4 vor Bagreagen: ee naqut 


TUE Sot nese 7 Sai ith tala bhve se { , 7 af G 
a7 : a a 


93 


Molz et al. (1968), modelling the effect of a single 
plant on water transfer, introduced the concept of potential 
soil-moisture availability which they defined as a measure 
of the capacity of a soil to transmit water to a root site. 
They developed a differential equation describing radial 
flow of soil water to a single vertical sink or root. The 
soil moisture flux was determined from this equation using 
Euler's method. 

In contrast to the microscopic approach - described 
above - of simulating the influence of plants on soil water 
transfer (Molz et al. 1968), there is the macroscopic 
approach (Molz and Remson 1970, 1971). In this case water 
transfer and extraction from an entire root zone is 
considered. Molz and Remson (1970, 1971) introduced an 
extraction term to the one-dimensional form of Richards’ 
equation and used the Douglas-Jones predictor-corrector 


finite difference method to solve the modified equation. 


Saturated and unsaturated flow in an integrated 
subsurface flow system 

Luthin and Day (1955) were among the first to model 
unsaturated and saturated flow systems as a single entity. 
They induced lateral flow through a sand-filled tank by 
maintaining a small difference in head between constant 
reservoirs. Flow nets were obtained by means of Laplace's 
equation, in which K=K(x,y), and which was solved by a 


relaxation method of numerical analysis. 


eloni=. 5s io 3902%9 gd3 paid bgabon },(8o8r): | 


Darsreyog to saaan03.ed2 besunongn’ panees 9 
ssdesen e as-fseiiep vad J ao priv anak talaiehe 


.93% 


aaseload do 2agem’ vd’ beittsads. ‘= as, 


fseroey ooiditewsan, noisasips isisheisesse 


sit doo) No Snte Lepissey-efotre 6 og¢tegee 


fav ROL ISBuUNS 2 ras cI. -perniniscan  2eFr nave 


[sbom’ O47 tear? ea enone a5 fh (tere ‘Yea 


qd «Ane? Caesltieonae? 6) dpuogss, wal fre 


+ +ayew Plechess’ 64: Loe e386 


Vv 

= 
c 
aX 


moet: noktdersxg, & 
if BasubSesseni (te | Seip Bie \Seeih 6 
eHosota 250 gee} Senos ae nth ie adinonaso9 a 
992599750801 9 38 esngl «ch guians| ad Soe , 
norseuos « sid igen ohz aVioe Be Sonuea 2 r 
Senta 


batez pasni 7S Gr WoL | inet 


7] 


io 
is 
a 
bik 
a 

p 
o-! 


; , 


\ 
' 
+ 
4 


Longe -€\ 2608 masega HeL2 badesn dine 
pens ¥ 


: { ; P yarns 2 = 
inedengs Assnded Deed at sees Lome 
7 7 7 


| 1" ie sd 7 
& Yo Bevioe eae igh Bnd ona 


os 


2 


The integrated approach to unsaturated-saturated flow 
systems was also used by Bouwer and Little (1959). In this 
case solutions to the two-dimensional steady flow problem 
were obtained by use of an electrical resistance network for 
relaxing the system. 

Freeze (1967, 1969) used a mathematical model to 
examine the mechanisms which control the relations between 
the unsaturated flow processes of infiltration and 
evaporation and the saturated processes of recharge and 
discharge. One-dimensional vertical, unsteady flow problems 
involving saturated-unsaturated systems were solved by means 
of a numerical difference technique. 

The next important advance in the continuum approach 
occurred when Rubin (1968) solved the Richards equation for 
two-dimensional transient flow in rectangular unsaturated or 
partly unsaturated soil slabs. The equation was solved by 
means of alternating-directions, implicit difference 
methods. 

Other workers who have featured prominently in the 
development of the integrated saturated-unsaturated flow 
System concept include Taylor and Luthin (1969), Hornberger 
et al. (1969), and Verma and Brutsaert (1970). 

The most comprehensive treatment of the unsaturated and 
Saturated flow systems as a continuum is undoubtedly that of 
Freeze (1971). He devised a three-dimensional model that 
simulated saturated-unsaturated transient flow in small 


nonhomogeneous, anisotropic basins. The model can be 
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collapsed to simulate two-dimensional flow, and can also 
accommdate hysteresis. Solution of the Richards equation was 
obtained using an implicit iterative method called the line 
succesSive overrelaxation (LSOR) technique. Freeze's model 
was discussed in some detail in the preceding section. 

Details of many of the numerical methods referred to in 
the foregoing discussion can be found in Remson et al. 


1971). 
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III. STRATEGY FOR MODEL DESIGN 
As stated previously, the primary opaecei ee of this study 
are: a) to develop a subsurface flow model in which 
vegetation appears as an integral part, and b) to use the 
model to simulate the effects of a variety of vegetation 
patterns on soil water and streamflow. 

The simplest approach is to consider the subsurface 
system as the core, and other process components, such as 
infiltration, evapotranspiration and seepage, as inputs and 
Outputs of the hydrologic simulation model. Inflows are 
regarded as gains or additions to, and outflows represent 
losses from, the system. A schematic diagram of the 
conceptual model is shown in Fig. 9. This simple model must 
relate to the physical watershed system. 

Given the complexity of soils and geology of most 
watersheds, a two-dimensional profile provides the best 
€onfiguration for simulating subsurface flow (Fig. 10); 
Two-dimensionality implies: that a given section has unit 
width, say 1m. A section so defined for simulating 
subsurface flow is assumed to supply generated streamflow to 
the adjacent 1 m of stream channel running perpendicular to 
the section. If the section is expanded uniformly in the 
third dimension, then the total outflow from this uniform 
"slice" is obtained by multiplying the outflow from the unit 
section by the length of channel adjacent to the slice. A 
simple example is a long uniform valley formed by 


escarpments (Fig. 11). This geomorphic pattern iS Often 
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encountered in the Front Ranges of the Rocky Mountains. 

The process of extending a two-dimensional problem to 
an entire three-dimensional watershed is fraught with 
Brcciculties. Conceptually, though, it is simple and 
involves four main steps: 

a) divide the watershed into several slices or segments 
based on geology, soils, vegetation, 
precipitation-elevation relationships, and other 
factors: 

b) carry out subsurface flow simulations using the 
two-dimensional, unit width representation (section), 
and generate outflow; 

c) use results from b) to determine the outflow to stream 
channel from each entire segment; 

d) route generated channel flows through the watershed 
using a standard routing procedure. 

Part a) can be readily accomplished if suitable data 
and detailed maps of the watershed to be simulated are 
available. Part b), the successful completion of which forms 
the main thrust of this thesis, is the most difficult step 
to resolve. It entails a statement of assumptions concerning 
the physical situation; the derivation and solution of the 
equations of flow through porous media; and the development 
of a mechanism for simulating outflow from the seepage face 
to the stream channel. 

The first necessary assumptions relate to defining the 


limits of the subsurface flow system. A datum is selected to 
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coincide with a lower impermeable boundary or basement (the 
X-ax1S AB in Fig. 10) across which no flow takes place. The 
position of this boundary is determined, preferably, from 
geologic maps. Vertical impermeable boundaries BC and DA 
(Fig. 10) define the upslope and downslope limits of the 
subsurface system. The presence of a no-flow boundary at BC 
implies that the phreatic divide coincides with the 
topographic divide. If it is further assumed that water does 
not move transversely (in the y-direction), then it is 
evident that no positive or negative leakage can occur into 
or out of the segment. All water leaving the basin does so 
via the stream channel or the atmosphere. Water is free to 
move across the upper boundary (ground surface CD in Fig. 
10) in either direction. It can also be transferred from the 
System interior by means of plants, pumping wells, or tile 
drains. 

Plants, especially trees, consume large quantities of 
water: often the volume of water that leaves a watershed as 
evapotranspiration exceeds the total volume of streamflow. 
During transpiration water does not cross the upper boundary 
directly, but is 'piped' up the stems and lost to the 
atmosphere through the leaves. For simulation purposes, it 
is useful to draw an analogy between transpiring trees and 
pumping wells. In the latter case, water is withdrawn from 
the saturated zone and discharged above ground; the well is 
considered to be a sink in the system. Analytical methods 


are available to determine the effects of pumping wells or 
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Sinks on the groundwater reservoir. 

Plants extract water, usually from the unsaturated 
zone, by means of an extensive root system. Research with 
forest trees (Kramer, 1969) indicates that water may be 
absorbed by trees at distances of several metres from their 
stems. Veihmeyer and Hendrickson (1938, cited by Kramer, 
1969) found that water 6 or 7 m from fruit trees is 
absorbed. The rooting depth varies with plant species. In 
most cases, the root zone is considered to be the upper 1 or 
2m of soil. However, some phreatophytes such as salt cedar 
(Tamarix chinensis Lour.) are capable of extracting water 
from depths of 6 m (Horton and Campbell, 1974), while the 
roots of mesquite (Prosopis juliflora) can penetrate as deep 
as 18 m into alluvium (Kearney and Peebles, 1951; cited by 
Horton and Campbell, 1974). 

In order to simulate water withdrawal by trees, a 
Series of sinks is introduced into the model near its upper 
boundary. The location of a sink in a finite element mesh 
Borcvesponds to a tree “rooting ‘depth ‘of’ 1 “to°2 m in’ the 
physical system. Each sink has a rate-of-withdrawal value 
assigned to it and represents the influence of a number of 
trees on subsurface water. When trees are harvested, 
simulation is achieved by deleting a suitable number of 
sinks from appropriate locations. Conversely, when regrowth 
occurs the sinks are restored. However, the 
rate-of-withdrawal value assigned to sinks representing the 


young trees will be less than that assigned to sinks 
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representing mature trees before harvest. This method allows 
for considerable flexibility since the differences in water 
demand due to species, age, and forest management practices 
are Simulated. The basic assumption here is that trees can 
be simulated using sinks in the manner described. The main 
limitation is the lack of information regarding consumptive 
use by different tree species at different ages. 

Although it is recognized that great differences in 
properties may exist within the same soil type or geologic 
Stratum, it is extremely difficult to map these spatial 
variations. Consequently, for simulation purposes, it is 
assumed that each soil type or geologic stratum is 
homogeneous. 

The remaining assumptions concerning the model are 
related to its mathematical development. In the saturated 
region, flow is laminar and Darcian; inertia forces, 
velocity heads, temperature gradients, osmotic gradients, 
chemical concentration gradients are all assumed to be 
negligible. Water, the soils and geologic formations are 
assumed to be incompressible. In the unsaturated zone, it is 
assumed that the soils are non-swelling and that the air 
phase is continuous and always in connection with constant 
external atmospheric pressure (Freeze, 1971). The 
implications of these assumptions are alluded to in the last 
Section (IIE) of the preceding chapter. 

The general equation of flow through porous media 


(equation 2.8) used in simulating water flow through a 
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hillside section is usually solved by application of 
numerical techniques. Such techniques involve a) 
discretization of a continuous domain, and b) generation of 
a set of linear algebraic equations. Discretization is the 
procedure whereby a continuous domain D is divided into a 
number of subareas for the purpose of solving ordinary or 
partial differential equations. Approximations to a 
continuous solution may be defined at isolated points by 
finite differences or, alternatively, defined over the 
entire domain by application of the finite element method. 
The solution in this case is the pressure potential at each 
point defined in the domain. 

The task of obtaining solutions for the parabolic 
partial differential equation of flow through porous media 
requires consideration of both initial and boundary 
conditions. The initial condition is the value of the 
function for the posed problem at starting time, to. The 
boundary conditions are either the value of the function, 
its normal derivative, or a linear combination of the 
function and its normal derivative on the boundary (Pinder 
and Gray, 1977). When the value of the function, such as 
constant pressure potential, is known, the boundary 
condition is referred to as a Dirichlet type boundary 
condition. If the normal derivative of the function is 
given, then the boundary condition is known as a Neumann 
type. Thus a vertical, impermeable or no-flow boundary 


forming a watershed divide exemplifies a Neumann type 
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boundary condition since dw/dx = 0, i.e. the hydraulic 

gradient in the x-direction is zero, so the flux across the 

boundary is also zero. 

Since the equation of flow for this problem cannot be 
solved easily by analytical means, both the finite element 
and finite difference methods will be employed to develop 
and solve the necessary approximating equations during 
Simulation. The finite element approximation will be applied 
to the spatial aspect of the problem, while the time 
derivative will be replaced by an appropriate finite 
difference representation. 

The finite element method has several advantages over 
finite difference techniques when applied to a continuum 
(Desai and Abel, 1972; Remson et al, 1971; and Rodarte, 
1978): 

1. It is not necessary to develop special expressions to 
define the boundary conditions. Close agreement with 
complex real boundaries are obtained. 

2. The dimensions of elements can be fixed arbitrarily. 
Thus, large elements can be used where a function 
changes gradually. When a function changes rapidly, such 
as in the near-surface regions of a simulated watershed, 
small elements are used. Irregular and complex 
configurations can be analyzed or simulated quite 
eaSily. 

3. Anisotropy and hetereogeneity can be routinely 


considered in the calculations. 
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4, It 1S possible to set up equations that produce 
positive-definite matrices which can be reduced to band 
form and solved with a minimum of storage and 
computation time. 

Application of the finite element and finite difference 
methods to the governing equations should yield solutions 
for simulated time intervals of, say, one day (since we are 
primarily interested in snowmelt and rainfall runoff events, 
and seasonal recession flows). The field of pressure 
potential which constitutes the solution at time, t, conveys 
information about the direction of flow within the profile, 
and the position of the water table. An additional 
algorithim is required to provide values of outflow 
generated at the seepage face. For this purpose an iterative 
scheme, developed by Neuman (1973) and described in the next 
chapter, will be used. 

The definition and solution of a simulated watershed 
section problem requires information regarding the geometry 
of the watershed, the hydraulic properties of the media, 
boundary fluxes, sinks, and streamflow data against which 
generated outflows can be compared. Thus, input data for the 
model will include some or all of the following: 

a) finite element mesh geometry of the entire vertical 
section, 

b) output from a snowmelt model to provide positive upper 
boundary fluxes during spring snowmelt periods, 


c) hydraulic properties of soil and geologic strata, 
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oye) amount,..intensity,.and duration. of rainfall, 

e) climate data required to compute evaporation, 

f) values of consumptive use by trees, based on species, 
age and perhaps other factors such as diameter at breast 
height (dbh), crown cover density, or height; 
consumptive use values are assigned to sinks located at 
near-surface points in the finite element mesh, 

The foregoing information provides both the initial and the 

boundary conditions of the problem. 

Output from the model will consist of daily values of: 

a) pressure potential (W) at points defined throughout the 
entire domain, 

b) total potential (¢), 

¢) volumetric water content (06), 

d) water table position, 


e) outflow from seepage face (generated streamflow). 
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IV. MATHEMATICAL DEVELOPMENT OF SUBSURFACE FLOW MODEL 


A. Governing Equations 
The full general equation for three-dimensional 
transient, saturated-unsaturated flow through porous media 


can be written (Freeze, 1971) as 


@_[go7(W)k,,(F, ae [go?(w)kyy(F, 3H) 


OX u Ox Oy u 


Hor gpo (yw) ks 2 try wie sasiMd (451) 
OZ 


OZ OM 


=[o(W)O(F,W){a' (F)te(F,W)8'(F)}+p(V)C(F,W) ] 
e(F, 7) se 


where 


ny), Zz coordinate directions, L; 

pressure potential, L; 

acceleration due to gravity, LT? 

viscosity of water, FTL™’; 

density of ‘water, FT*L *; 

intrinsic, or specific permeability, L’; 

porosity, dimensionless; 

volumetric moisture content, dimensionless; 

geologic formation or soil strata, dimensionless; 

=00/d~ specific moisture capacity, L™'; 

"=apg vertical compressibility of formation, L™' 
coefficient of vertical formation compressibility, 
ja 

=Bpg compressibility of water, L™' 

coefficient of water compressibility, L’F™' 

times Ts 


Fk aADTWA FTDENO EX 


TRMD 


If the foregoing expression is written in vector notation it 


Simplifies (Pinder and Gray, 1977) to 
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V°[Kk,.* (Vytk) ]=Cawt+Os, av (4.2) 
ote ot . 
where 
K=pgk,,;/u saturated hydraulic conductivity tensor, LT-'; 
Kew relative permeability of water, dimensionless: 
k unit vector in the vertical direction: 
Ss, specific storage, L™'. 


B. Finite Element Simulation 

The finite element method was developed in the 1950's 
to solve complex structural engineering problems. Since that 
time it has been successfully applied to problems in: soil 
and rock mechanics, heat conduction (Desai and Abel, 1972), 
surface and subsurface hydrology (Pinder and Gray, 1977), 
and in simulating water uptake by plants (Feddes et al, 
1976). It has also been used to describe the geometry and 
Dosition of cells in fossil and living plants (Niklas, 
977) 

The following development of A Gaaiie element model for 
simulating unsteady, saturated-unsaturated flow of a single 
fluid through porous media is taken from Pinder and Gray 
(1977), and from Neuman (1973). The basis for the finite 
element model is Galerkin's method of weighted residuals. It 
is essentially a scheme for solving differential equations 
using integral approximations. We first consider the time 
independent equation 

Lu = f in B (4.3) 


where u is an unknown function, f is a known functron, Uo1s 
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the operator, and B is a bounded domain. Let u be 
approximated by a function i(x) which consists of a linear 
combination of suitable functions and satisfies appropriate 
boundary conditions of the boundary value problem. For 


example, we might try 
Pe M 
O(x)=,(x)+Za ,¢; (x) (4.4) 
j=2 
where ¢,(x) is chosen to satisfy the essential boundary 
condition and the coordinate or basis functions ¢,(x) (also 
called bases) satisfy the corresponding homogeneous boundary 
Sendition., Usually ¢;(x) ts not written explicitly, but is 
incorporated in the series 
. M 
G(x) =Za ;; (x) (4.5) 
j=1 
Let a residual R be defined by 
R(¢x)=Lu-t (x) (4.6) 
Inserting the trial function, we have 
Ltt (4.7) 
R(x)=L[Za ;¢; (x) ]-£(x) . 
JS) 
The residual (or error) function R(x) will be zero for the 


exact solution, but not for the approximation. In the method 


of weighted residuals, the residual is minimized by settings 
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the inner product of the residual function and a weighting 


function w,;, equal to zero: 


f, R(x) ew, dx=0 Wai 29 <M) (4.8) 


This equation can more conveniently be written as 


<R(x),w,>=0 ira tn MY (4.9) 


There are several weighted residual techniques which 
can be used to obtain the solution to a set of differential 
equations. They include the collocation, subdomain, least 
squares, and Galerkin's, methods. The method selected 
governs the choice of the weighting function, w;. In the 
Galerkin method, which was selected for this study, the 
coordinate functions ¢,(x) also serve as the weighting 
functions w,(x). When the coordinate functions are inserted, 


the previous equation becomes 


<R(x),¢;,>=0 Cea e2e eee) (4.10) 


In order for the residual to be minimized, R(x) must be 
Orthogonal to all the functions ¢,(x). 

The principle of the finite element method is that the 
coordinate functions, which are usually polynomials, are 
piecewise continuous over subdomains called finite elements. 


Nodes are located along the boundaries of each element and 


gticsipiow 5 ben nol gongs: dinate ‘an 


“apts 49 1 


-®) 


S) (My vino Bebe 7 Qariy 


(M4 Sale ear? ‘ One ae a 
Aaa) 
! 
2a,plotoss J subiege bsiipiewlerervem aie 
Dib ta seas OF Pelaoion “ae4), mises 
‘onol pc amobdue (nebgeselie>oneteeiriaat, @ 


2.32 DonISsh; sar ichabat t 1% HEsretad | aT: 


Gl w \fogtonbarpalgnpe wba erions 


- a 


? 
, Se 
ei: . ybUse(sind ares base si ae Ri ; 
Snr sibs: SA? se a¥% ‘ae \ sath cone snot; 
r ae | 
~ ee aie ¢ +, ec 7 . i 

920) 845 zadiyonee Seen eeoee ei? ech Fe 

ei ae out 12 

eamiode debs avery 

ad yey? ] baci b a 

/ . 


(at .3) ) aicdgeetiny S ) idee pa lies 
— Ne Se 
oa) An 
ae 
| as “ 
ht afckvonsi aid Pie-® 


ad ‘taum cx) 8) cesta ns 63  lsubigat 


add -té6dj-8h boise jnpmiets etinia $a, tp aan 
ste .elaiiio agheg tecaa hide Hatt - ‘ 


‘2icemele ot2nit baligo ey 


110 


each coordinate function is identified with a specific node 
i.e. $,(x) has a value of unity at node i, but goes to zero 
at all other nodes. The elements may be one-, two-, or 
three-dimensional; linear, quadratic or cubic in space 
and/or time; rectangular, triangular or tetrahedral in 
shape. They may have either straight or curved sides. These 
examples are only a few from a large number of possible 
elements. 

The preferred element for this model is the linear 
triangular element because it is well suited for use on 
irregular boundaries and concentrates coordinate functions 
in those regions of the domain where a rapidly varying 
solution is anticipated. Coordinate functions may be 
formulated in terms of either global or local coordinates. 
For triangular finite elements the local (or area) 
coordinate system is used because element integrations are 
thereby greatly facilitated. The development of coordinate 
fuctions using area coordinates for triangular elements and 
integration of the resulting expressions are described in 
Appendix A. 

The general expression for coordinate functions for 
linear triangular elements is 

gd, =a, tb; x+cj\y Carly) 
where a,, b,;, and c, are constants identified with the ith 


node and the ith coordinate function. When they are 
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expressed in terms of triangular area coordinates, L,, the 


coordinate functions become 


Gy XY okey) t(Yj-y.) xt (xox) yd 2A. a ear-%) 
@ =L =((x.yirxiryn) + lyenyi )x+(xi-x, )yJ/2a. (4,12B) 


Oem mux; ¥;-X7y.) ly poy, )xt0xyox, yl 2A ace) 


where A,is the element area. 
The foregoing procedures will now be applied to the 
governing equation, 4.2. Let the pressure potential, Ww, be 


approximated by 
N 
w~W=zP , (t)¢; (x,y) (04.94:3)) 
Jel 
The two-dimensional form of equation 4.2 is 


r 


Vay (Rk Vayw)t¥yy (Rk ok) -[C+8S, ]ay=0 (aves) 
e Che 


and Galerkin's method requires that. the coefficients P be 


determined such that 
<LY,¢,>=0 (Bias eyye eet ND (4,15) 
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=A os ee v K 
<Vxy° (Kk, Dh (t)¢; (x,y) ,6)>+<Vy,° (Kk 


N 
~<[C+O0S,]0(ZP, (t)o,(x,y)) ,¢6;>=0 (4.16) 
e J=1 Ot 


Application of Green's theorem eliminates the second 


derivative, and reduces the governing equations to the form 


<K' (0) Vx yP 505 ,VxyGi>t<[Ct+OS,]¢,dP,o;> 
e dat 


=f FeR(b) LV. d+k] 9 .ds (4.17) 
where K'=Kk,, and fi is an outwardly directed unit normal 
vector on the surface, I. 


This equation may also be expressed in matrix form as 


[A]{P}+[B]{dP}={F} (4.18) 
dt 


where the elements of matrices [A], [B], and {F} are 


aj r=<K' i) Vio Ne) O poy 


b, ;=<(C+0S.]¢,,¢i>, 
e 


Ei, =J A-K' (0) -[V.,0+k]¢,ds 


In general, for the subsurface flow problem, [A] is 
Symmetric and [B] is a diagonal matrix. 
The earlier assumption that the compressibilities of 


water, soil, and geologic formations are equal to zero 
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causes the specific storage, S,, also to become zero. Thus 


b;,; reduces to 
b; ;<[C]¢;,¢;> (4.19) 


Considering first the steady state problem, equation 4.18 


reduces to 
[A]{P}={F} C4720) 


When triangular finite elements are used, the information 
obtained from integration over each element is stored ina 


3x3 element coefficient matrix as 
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The components of these matrices are solved using the 
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following relations 


<09; , 99 j>+<06; ,00;>=1_ (ly j-yn) (ye-yi) + (xn-x |) (x5 -x_)) C475 2)25) 
Ox ax dy doy 4A, 


and 


<0¢ i 1 9b, >+<09; 19 
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ieee (Cys aye) etx x) (4.23) 
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The next step is to assemble the element information 
into a global matrix, [A], by summing, for a given node, the 
contributions to that node from each element coefficient 
matrix. If the transient problem is considered, the global 
matrix [B] must also be assembled. This is done in much the 
Same way as was matrix [A]. However, the element coefficient 


matrix in this case has the form 


<61,01> <b1,62> <$1,93? 
<h2,91> <h2,62> <%2,637 (4.24) 
<G3,917 <b3,92> <b3,937 
The integrals are solved by means of the expressions 
<9, 19) >=A_/6 (l=i,j,k) (4.25) 


and 


<9) 1om>=0 (l=i,j,k), m=i,j,k; 1ém (4.26) 
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C. Application Of Finite Difference Approximations To Time 
Domain 
There are several methods available for replacing the 
time derivative in equation 4.18 with a finite difference 
approximation. The most general is the weighted average 


approximation 


bat) Cai Be Sh eiser) iP} i+ hela) Gf P hNaoe fe p) 
At 


=e{F}t* 4t+(1-e) {F}* (O0<e<1) een. 


where t is the time, and At is the time step. 
When e=0, an explicit finite difference scheme results. This 
method is the simplest to solve but is only conditionally 


Stable. The explicit scheme is 


GAS ekB) ) i PL eT EB ip} ot={ hh (4.28) 
At At 
An unconditionally stable, first-order correct scheme is 


obtained when e=1, and equation 4.27 becomes 
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This scheme is usually referred to as the implicit or 


backward finite difference method. 


A time-centred scheme results when e=1/2. Commonly 


referred to as the Crank-Nicolson implicit method, it~ is 
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second-order accurate and usually neutrally stable. It can 
be written as 


LAvCER STP} ©) [BIC {Pp} ** A= {Pp} *) 


1 +1 
Z At 


Sete ee poe) (4.30) 
2 


Because of the dependence of both K' and #6 on pressure 
potential, equation 4.18 is highly nonlinear, and the set of 
nonlinear algebraic equations it represents must be solved 
uSing an iterative procedure. A solution may be obtained by 
combining the finite element method with either the 
time-centred or the backward difference scheme. In this 
case, a backward difference scheme for time-centred 
coefficient matrices is used. (Pinder and Gray, aaioa7 > 
Neuman, 1973). It is effectively an under-relaxation 
technique and has the advantage of damping oscillations that 
frequently occur when solving highly nonlinear systems. This 
iterative scheme also overcomes the problem of obtaining 
solutions that arises when the specific storage, S., 1S 
zero. 

The fully implicit backward difference formula for this 


problem is 
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where t represents the time, and At is the time step. 


Equation 4.31 may be expressed as 
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where 


Mod mac? Ate GAlt i GR Aet+1 [Bats ile At 
At 


ee hs Ate Sire oh ph 
Bt 


and m is the iteration number. 

Several steps are required to complete the iterative 
procedure which provides the solution to this equation. 

1. The initial conditions for the problem are described at 
each node. They include initial volumetric water content 
6=0(W), and initial hydraulic conductivity K'=K'(y). 
Thus initial pressure potential, w (and {P}*), is also 
prescribed. 

2. An initial estimate of pressure potential at the end of 
the-first timesstep, -{P}!*“', is obtained by ‘solving «the 


expression 


[GE] {ip} *+4t={p} t+{Fr}* €4°. 3:3:) 


rather than solving equation 4.32. 


3. This estimate is used together with the relation 
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{P} tt '/2Ata({p}t+{p}t+ At) /2 (4.34) 


to determine {P}'*'/24t 

The values of {P}**'’?4* govern the values of K', C, and 
6 which appear in coefficients a;, and b,,; (equation 
4.18). New values of K, C, and 6 are determined from 
{P}t*+'’24t and the resulting updated coefficient 
Matrices are used in equation 4.32 to obtain an improved 
valuesforoe{P}**A*. 

The iterative procedure (steps 3 and 4) is repeated 


until the difference between successive iterations is 


within satisfactory limits i.e. 


er Ate tc) eB, t+ Atm <5 (eile zee gee NN) (4.35) 


where E is the specified difference criterion. 
When this condition is met, {P}* is updated to {P}*t*“', 
and the iterative cycle is initiated once again for the 
next time step, solving first for {P}**?4*-". 
For constant-sized time steps beyond the initial time 
step, the following relation is used to obtain a first 
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The matrices generated by using the finite element 
method are usually sparse and banded, a condition which 
permits considerable saving in core-storage and 
computational effort (Pinder and Gray, 1977). There are 
several methods available for reducing matrices to a form 
Suitable for extracting the solution, {P}. The method 
selected for this problem entails column by column 
decomposition of the matrices, and is referred to as the 
Cholesky method. The theory, procedure, and the algorithms 
necessary to implement the Cholesky method are given in Elwi 


and Murray (1977). 


D. Treatment Of Boundary Conditions, Sinks, And Seepage 
Faces 

The boundary conditions are handled through the line 
integral in {F}. This integral is non-zero only along the 
domain boundary, I, where it represents the component of 
flow normal to [. For a no-flow boundary condition, {F} is 
set to zero at each of the relevant boundary nodes. If {P} 
is known for a given node (Dirichlet-type boundary 
condition) then the equation for that node is condensed or 
partitioned from the matrices. 

It appears that no provision for sinks has been made in 
equation 4.18. Strictly speaking a term, {OO} EOrasInks 
should appear on the right hand side of the equation. 
However, for computation purposes {Q} has, in effect, been 


absorbed by {F} such that 
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{F}'={F}+{Q} (a3 7) 


Since {F} is zero for all interior nodes, and {Q} has 
non-zero values only at specified interior nodes, it is 
feasible to use {F} to define both {Q} and {F} without 
changing the definitions and assumptions related to the 
vector, {F}. 

Consideration must also be given to proper simulation 
of outflow from seepage faces during iteration. For this 
eee the approach devised by Neuman (1973) will be used. 
This method predicts the location of the seepage face at 
time t"*', given its position at t". It works as follows: at 
time, t", the seepage face is defined and the pressure 
potential, P, is set to zero at all points along the seepage 
face. At points where P<0, the inward-directed normal flux, 
Q, is prescribed as zero. The solution to the governing 
equations should then yield Q<0 (outflow) when P=0, and P<0 
when Q=0. If, however, results show Q>0 at a point where 
P=0, then Q is set equal to zero, and the point treated as a 
no-flow portion of the boundary. Similarly, if P>0O where 
Q=0, the corresponding value of P is set to zero and the 
point treated as a constant potential point during the next 
iteration. errr is continued until results are 


compatible. 
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V. SIMULATIONS: RESULTS AND DISCUSSION 
A. Unsaturated Flow 


Model validation 

The credibility of a physically based mathematical 
model should be established before it is used to simulate 
the complex, dynamic systems of the real world - systems in 
which not one, but several processes may be active at the 
Same time. One has first to ask himself the question: "Is 
the model responding correctly, in a qualitative way, to the 
boundary and initial conditions imposed upon it ?" Applied 
to a specific hydrologic process, the question might be 
re~phrased as: "Does the model respond with reduced water 
content, lower potentials, and steep potential gradients, 
near the ground surface during times of high evaporative 
demand ?" Alternatively, one might pose the question: "Can 
the model show an advancing wetting front, and a decay 
muInction fOr Infiltration srate, during) Infiltration ?°°" Such 
relationships were established Vong agc (Horton 1933, for 
example) through extensive field observations and laboratory 
experiments; they are described in many soils and hydrology 
text books. 

Once it has been determined that the physically based 
mathematical model can, indeed, simulate isolated processes, 
then it is feasible to use the model to examine such 


processes within the context of an integrated system, such 
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as a watershed or a hillslope. 

It is the intent, in the following simulations, to have 
the subsurface flow problem evolve from a very simple case 
Gieinfiltraticn; evaporation, or transpiration for a single, 
unsaturated medium, to more complicated situations in which 
Saturated-unsaturated flow systems, several media, and the 


seepage face are considered. 


One-dimensional, vertical infiltration and evaporation 

The first simulation concerns one-dimensional flow 
through unsaturated clay soil. This problem was considered 
first because it is very simple. The volume of soil involved 
is small enough that the usual assumptions, that the soil is 
homogeneous and has uniform properties, are probably correct 
in this case. Furthermore, a solution to the problem has 
been obtained by other investigators (Beven, 1975; Philip, 
1957). Philip solved the problem using an analytical 
technique, and Beven arrived at essentially the same result 
by means of a finite element method. 

In each case, flow through Yolo Light Clay was 
simulated. This soil, which is an agricultural soil in 
California, consists of 23.8 percent sand, 45 percent silt, 
and 31.2 percent clay, with a pore space of 50 percent 
(Moore, 1939). The unsaturated hydraulic conductivity is 
related to pressure potential by the expression 
K=0.00463/(w?7+400.0), where yw is the pressure potential 


(Rubin, 1968). Thus, saturated K is 1.16 X 10°* cm/sec. 
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A mesh of unit width with no-flow boundaries on either 
side (Fig. 12a), suggested by Beven (1975), was constructed 
for the one-dimensional problem. Note that z-coordinates for 
all but two of the nodes are negative.This has no effect on 
the total potential at one node relative to the total 
potential at other nodes. The only difference between 
results from a positive coordinate system, referenced from 
the bottom of the soil column, and the results from this 
System, is that the total potential at all nodes in the 
negative system is 10 cm higher than the potential at the 
corresponding nodes in the positive system. Initial pressure 
potentials were set everywhere equal to -660 cm (6=0.240), 
except at the surface nodes where they were maintained at 
Saturation (W=0 cm). These initial and boundary conditions 
correspond to those of Beven and Philip. 

The results produced by SUBFEM at the end of 10,000 
sec, using a time step of 500 sec. (Fig. 12b) compare 
favorably with Beven's (1975) finite element solution and 
Philip's (1957) analytical solution. In all three cases 
water penetrated to a depth of about 8 cm. Some differences 
are inevitable since the functional relationship between 6 
and Ww, and between K and yw, used in SUBFEM are different 
from those used in Beven's and Philip's simulations. 
However, the retention curves do have certain data points in 
common. Convergence occurred fairly rapidly in the 


one-dimensional simulation - to a difference in potential of 


Me01 cm after 4 to 8 iterations. 
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b) 


VOLUMETRIC WATER CONTENT (@) 


0.24 | 0.496 


t= 10s sec 


Analytical solution (Philip, 1957) 
—x—x-= Finite element solution (Beven, 1975) 


—oO——O— SUBFEM solution 


Finite element discretization for one-dimensional 
infiltration into Yolo Light Clay (After Beven, 1975). 
Comparison of SUBFEM solution with those obtained 
by Philip (1957) and Beven (1975). 
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The isopotential lines for the soil column, after 
10,000 sec, are shown in Fig. 13a. The highest potential (0 
cm) occurs at the upper boundary and the lowest (-660 cm) at 
a depth of 10 cm. This confirms that vertical infiltration 
mS taking place, i.e. in the direction of decreasing 
potential. The flow lines or streamlines run perpendicular 
to the isopotential lines 

Field capacity is defined as the water content of soil 
after gravity drainage is complete, usually at a tension of 
1/3 atmosphere, or -343 cm of water. After 10,000 sec of 
infiltration, the soil column is at field capacity, or 
wetter, down to a depth of about 6.7 cm. The lowest point at 
which field capacity is attained is in the region of 
steepest potential gradients - between depths of 6 and 8 cm 
(Fig. 13a). It is also in the region where the greatest 
changes in water content occur (Fig. 13b). 

The zone where water invades and advances into 
Originally dry soil is known as the wetting front zone. It 
is the zone in which the greatest water-moving forces, 
resulting from potential gradients, are produced (Hillel, 
1980). For this problem, it corresponds to depths 6 to 8 cm, 
where the isopotential lines are close together and the 
potential gradient is 199 cm/cm (Fig. made, eimieonGrast:, the 
isopotential lines between 0 and 6 cm depth are more widely 
Spaced, yielding a potential gradient of 34 cm/cm. 

In order to evaluate transient behaviour of the 


simulated flow, the position of the wetting front is plotted 
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Vertical exaggeration :0-5x horizontal scale 
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(a) TOTAL POTENTIAL (cm) (b) VOLUMETRIC WATER CONTENT 


One-dimensional, vertical infiltration - total 
potential (a) and volumetric water content (b) 


after 10,000 sec. 


Pagure 13. 
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over time (Fig. 14). Here, the wetting front is defined 

eet ani ly as the depth at which the water content is at 
field capacity. The curve shows a smooth progression of the 
wetting front from the surface at t=0 sec, to a depth of 
nearly 7 cm after 10,000 sec. This is the "advancing wetting 
front” referred to earlier. 

If the flux across the upper boundary is plotted over 
time, the resulting curve can be described by a decay 
function (Fig. 14). This is the infiltration capacity curve 
Por the soil, and theoretically should become horizontal 
(constant) when the infiltration rate equals the saturated 
Pearaulic conductivity: ct the soil, (1.16 X 10> *ecm/sec). 

A transient, one-dimensional evaporation Simulation was 
also run. It differed from the infiltration run only with 
respect to the boundary conditions at the surface nodes. A 
negative flux corresponding to an evaporation rate of 5 
mm/day was specified for these nodes, instead of the 
‘constant pressure potential (w=0) condition that applied 
during infiltration. The choice of 5 mm/day is reasonable as 
it is a rate commonly attained on warm, windy days in the 
mountains and foothills of Alberta. The initial conditions 
remained the same as for the infiltration problem. 

At the beginning of the simulation, the water content 
ranged from 0.496 (saturation) at the surface to 0.250 at 
0.5 cm depth. For the same depth increment, the pressure 
ranged from 0 to -660 cm. This is equivalent to a potential 


gradient of 1320 cm/cm in the downward direction. 
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Consequently, when the simulation proceeds, there should be 
some evidence that drainage is taking place. 

After 10,000 sec, the combined influence of evaporation 
and drainage had produced a potential gradient in the upper 
1/2 cm of soil, of about 60 cm/cm. downward (Fig. 15a). This 
is considerably less than the gradient in the same depth 
increment at the start of the simulation. The gradient's 
downward direction indicates that drainage is the dominant 
process in this case. 

Both total potential (Fig. 15a) and water content 
(Fig.15b) show little change from initial conditions at 
depths below 4 or 5 cm. This is consistent with field 
observations which show that water draining from the 
surface, wets progressively deeper layers with time. 

At the beginning of the simulation, water content was 
uniform throughout the soil column, except in the surface 
increment. After 10,000 sec, the water content in this 
increment approached that of the rest of the column. For the 
soil column as a whole, the water content tends to decrease 
einightly with depth (Fig. 15b). It is not possible in this 


case to discern the effects due to evaporation. 


Two-dimensional, horizontal infiltration 

The two-dimensional, infiltration problem was selected 
because it, too, has been solved by others (Beven, 1975; 
Pein: 1968). In this case horizontal infiltration into a 
10-cm high block of Yolo Light Clay was simulated. The total 


hydraulic potential (@) at the infiltration face was 
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Evaporation : 5mm/day 


Vertical exag: x0-5 


Bs 10.0 L 
(a) TOTAL POTENTIAL (cm) — (b) VOLUMETRIC WATER CONTENT 


(cm?/cm?) 


Figure 15. One-dimensional evaporation ~- total potential (a) 
and water content (b) after 10,000 sec. 
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maintained at -13 cm, while the rest of the block was 
initially at w=-568 cm (Beven, 1975). The K vs w relation, 
and the water retention curve were the same as for the 
one-dimensional problem. A time period of approximately 6 
hours (22,000 sec) was simulated using Beven's (1975) finite 
element discretization (Fig. 16), and a time step of 1,000 
sec. 

Figure 17 shows the total potential field for this 
problem, obtained using Beven's finite element model (Fig. 
iva), Rubin's method (Fig. '17b)/ and SUBFEM (Fig. 17c). 
Rubin (1968) employed the alternating-direction, implicit 
difference technique to obtain his solution. The three 
solutions show quite good agreement - the equipotential 
lines for each are almost vertical, and indicate horizontal 
Peow from leftyto right. No difficulties related to 
convergence and numerical stability were encountered for 
this problem. 

A comparision of the block's volumetric water content 
(Fig. 18) determined by Rubin (1968) and SUBFEM indicates 
that after 6 hours the water content is about 3 to 8 percent 
higher for the SUBFEM simulation. The difference may be 
attributed to the 9 vs W relation used in SUBFEM, which is 
different from that used in Rubin's (1968) simulation. 
Another factor could be the slightly greater simulation 
period (6 hours 7 minutes) used in SUBFEM, compared with 6 
hours used by Rubin. In SUBFEM, the volumetric water content 


corresponding to the initial pressure potential (-568 cm) is 
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Figure 17. Two-dimensional, horizontal infiltration - total 


potential after 6 hr. 
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Figure 18. Two-dimensional, horizontal infiltration - water 
content [a) and b)] and total potential [c)] fields 


after 6 hr. 
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weeos, Dut in Rubin”™s simulation it was 0.238. Tt is also 
possible that a lower limit for the convergence criterion in 
SUBFEM would have produced solutions that were closer to 
those of Rubin. After 6 hours, water had advanced about 13 
emeento the block (Fig. 18). 

The soil is at field capacity at a horizontal distance 
@reapout il cm from the infiltration face (Fig. 18b). This 
gives us a reasonable indication of the position of the 
wetting front. The potential gradient between 10 and 12 cm 
is about 90 cm/cm, and between 0 and 10 cm, 29 cm/cm (Fig. 
fee). This relation, or pattern, is very similiar to 
conditions described for the vertical infiltration problem 
which showed steep potential gradients in the vicinity of 
the wetting front, and lower gradients near the infiltration 
face. 

The advancing wetting front curve for horizontal 
Mitaitration (Fig. 19) is very similiar to that .for) vertical 
infiltration CElg.e 14) ear tere iO, 000 sec a the wetting front 
had advanced 6.7 cm into the vertical column. The 
corresponding value for horizontal infiltration is 7.1 cm. 
The infiltration capacity curve for horizontal flow (Fig. 
19) resembles the one for vertical flow (Fig. 14) in that it 
portrays a decay function. 

The results obtained so far, for vertical and 
horizontal infiltration, suggest that SUBFEM can simulate, 
qualitatively at least, both isolated and interacting 


hydrologic processes such as infiltration and drainage. The 
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solutions are generally consistent with those obtained by 
other investigators and also with observations from field 
and laboratory experiments. The other researchers employed 
finite element, finite difference, and analytical techniques 
to achieve the same solution. On the basis of the foregoing 
sults, it would appear that, since SUBFEM can simulate in 
a limited way, real hydrologic conditions, its validity has 
been established in part. Thus we can proceed, with some 
measure of confidence, to use the model to try and solve 


more complicated problems. 


Evaporation and transpiration from a large box of soil 


Evaporation 

A run was performed to determine model response to 
evaporation simulation. In this case, evaporation was 
assumed to occur from Yolo Light Clay soil contained ina 
large, open box or container 30 m long, 1.8 m high) sand ot 
arbitrary width. The system might be considered as an 
extended lysimeter. A simple finite element mesh, consisting 
of 64 nodes and 90 elements, was used (Fig. 20). No-flow 
boundaries exist at x(z,t)=0, x(z,t)=3000, and z(x,t)=0 cm. 
Evaporation (negative flux) nodes are located at the 
surface, and an evaporation rate of 5 mm/day was specified. 
Initial conditions were set so that w(x,z,0)=-500 cm 


(6=0.260). The problem was run for a simulated time interval 


of 240 hr using a time step of 12 hr. 
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The results show (Fig. 21) that for this problem, 
effects of evaporation on the soil are limited to the upper 
60 or 70 cm. No changes were produced in the lower portion 
(heights 0 to 110 cm) of the profile. Pressure potential 
values after 10 days (Fig. 21a) ranged from -3330 cm at the 
surface to -500 cm at depth 70 cm below the surface. At the 
surface nodes, differences between pressure potentials over 
consecutive time steps tended to increase with time. 

Steep potential gradients, up to 40 cm/cm, exist in the 
upper 70 cm of soil (Fig. 21b). They arise in order to 
sustain the evaporation rate at 5 mm/day, compensating for 
the reduction in the hydraulic conductivity which occurs as 
the soil dries out. The gradients are directed upward which 
is consistent with the direction of water movement in 
response to evaporation. 

After 10 days of evaporation, the water status of the 
S0il has changed significantly (Fig. 21c). At the surface, 
the water content has been reduced by 10 percent. Below the 
Surface , reduction in water content decreases with depth 
for about 70 cm, below which point uniform moisture 
conditions (@=0.260) prevail. 

The boundary conditions, coupled with the fact that 
evaporation was not imposed on the surface corner nodes, 
have a pronounced effect on the shape of the isolines (Fig. 
21). The isolines bend sharply near the lateral boundaries 
and intersect the surface at right-angles. The total 


potential lines (Fig. 21b) indicate that water is drawn away 
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from these boundaries during evaporation. 

The isolines shown in Figure 21 display a symmetric 
pattern about the 1500 cm horizontal distance mark. 
Supplementary simulations suggest that, if the box is 
tilted, this symmetry disappears. Positions downslope from 
the 1500 cm mark then tend to be wetter, and positions 
upslope dryer, than the mid-slope positions. These 
Supplementary simulations also revealed that the proportion 
of the profile affected by evaporation increases with 
Saturated hydraulic conductivity. If the conductivity is 
great enough, water may be removed from the bottom of a 


profile, even during the first time step. 


Transpiration 

For the purpose of simulating transpiration, it was 
assumed that trees are planted in the box used in the 
previous simulation, and that the box is otherwise 
completely sealed to prevent soil evaporation. It was 
further assumed that the trees are of sufficient size and 
Number to transpire 5 mm of water per day. 

A few changes only are required to convert the 
simulation for the previous run from an evaporation to a 
transpiration problem. Transpiration is simulated by 
activating sink nodes at depth 60 cm below the surface, and 
assigning to them appropriate values that correspond to 


transpiration rates of 5 mm/day. Since no evaporation 15S 


Securring, fluxes at the surface nodes are set to zero. 
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The results indicate (Fig. 22) that transpiration 
affects a greater volume of soil than does evaporation CBG. 
21), under the same conditions. The volume affected extends 
down to about 110 cm below the surface. The lower 70 cm are 
not affected. After 10 days, the lowest total potential 
(~1447 cm) was attained at the sink nodes (trees), 60 cm 
below the surface. 

Because tree roots are completely surrounded by moist 
soil, they have access to water in every direction. In 
contrast, evaporation is controlled by surface conditions, 
and water has to be extracted from ever-increasing depths. 
As the surface soil dries out, it transmits water less 
readily, even if water is available at lower depths. 
Consequently, more water is removed from the surface layers. 
This is reflected in the low pressure potential (-3330 cm) 
values at the surface for the evaporation problem. The 
lowest pressure potential achieved for the transpiration 
problem was -1567 cm. Differences between the evaporation 
and transpiration problems are further reflected in the 
potential gradients. Evaporation produced gradients of up to 
40 cm/cm - iehey the average gradient (about 20 cm/cm) for 
the transpiration problem. 

The isopotential lines (Fig. 22a) are almost concentric 
about the row of sink nodes, which implies that water is 
drawn inwards from the surface, the lateral boundaries, and 
from the soil below. The potential gradients are slightly 


Steeper (22 cm/cm) below the sink nodes than they are above 
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TRANSPIRATION FROM A FLAT PROFILE (5mm /day) 


HEIGHT (cm) 


1000 2000 
horizontal distance (cm) 


(a) TOTAL POTENTIAL (cm) 


Vertical exag : x8 


1000 2000 3000 
horizontal distance (cm) Ses 
(cmifom?) 


(b) VOLUMETRIC WATER CONTENT 


Transpiration from a flat profile - total potential (a) 


Bicure 22% 
and water content (b) after 240 hr. 
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them (18cm/cm). 

Data for individual time steps indicate that 
transpiration did not affect the surface until after about 3 
days. 

After 10 days, the volumetric water content in the 
vicinity of the sink nodes had decreased by about 6 percent 
(Fig. 22b). The reduction decreases away from the sink nodes 
toward the lateral boundaries and to depth of about 110 cm 
below the surface. At these points there is little or no 
change from the initial conditions (6=0.260). There 1s also 
a reduction in water content above the sink nodes ranging 
from 6 percent at the sink nodes to about 0.4 percent at the 
surface. These data suggest, what is self-evident in the 
field, that during transpiration, trees remove water from 


the soil surrounding their roots. 


Evaporation plus transpiration 

The box described in the two previous experiments is 
opened up for this simulation to permit evaporation and 
transpiration to occur simultaneously. In every other 
respect, the three problems of evaporation, transpiration, 
and evaporation plus transpiration, are identical. 

Both the surface and the sink nodes are activated for 
the evaporation plus transpiration simulation. The 5 mm/day 
water withdrawal rate used previously is divided equally 
between evaporation and transpiration so that each occurs at 


the rate of 2.5 mm/day. 
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The plot of total potential for evapotranspiration 
(Fig. 23a) shows features that were evident when evaporation 
and transpiration were treated separately. For example, the 
lowest total potential (-1690 cm) occurs at the surface 
which was also the case in the evaporation simulation (Fig. 
21b). The effects of evapotranspiration extend to a depth of 
about 110 cm - the same depth as for the transpiration 
Simulation (Fig. 22a). 

The potential gradient produced by evapotranspiration 
(11 cm/cm) is only 1/4 of the gradient resulting from 
evaporation (Fig. 21b), and half the gradient resulting from 
transpiration (Fig. 22a). The isopotential lines indicate 
that flow is directed upward, and away from the lateral 
boundaries. 

The reduction in volumetric water content at the 
Surface is greater for evapotranspiration (Fig. 23b) than 
for transpiration (Fig. 22b), but less than for evaporation 
(Fig. 21c). The reduction decreases almost linearly with 


depth. 


Evapotranspiration from sloping profiles 

If it were possible to tilt the box to, say, an angle 
of 20°, and open the downslope end so that water is free to 
move out, we would have a very crude representation of a 
two-dimensional, unsaturated flow system for a forested 


hillslope. 


These added complexities can be accommodated within the 


context of the previous problem by simply changing the 
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COMPLETELY FORESTED 


Evaporation > 2-5mm/day 
Transpiration : 2:5mm/day 
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HEIGHT (cm) 


1000 2000" =, 3000 
horizontal distance (cm) 


(a) TOTAL POTENTIAL (cm) 


Vertical exag : x8 


HEIGHT (em) 


3000 


1000 2000 
horizontal distance (cm) eee 
(b) VOLUMETRIC WATER CONTENT (cm/crv ) 


Evapotranspiration from a flat profile - total potential 


Pieure 23, 
(a) and water content (b) after 240 hr. 
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z-coordinates and redefining the downslope (left) no-flow 
boundary as a negative flux (evaporation) boundary. 

Until now, for vertical and horizontal flow, the same 
coordinate system has been used to compute the solution and 
to plot the isolines. However, for sloping profiles the 
procedure is \different. The original coordinate system, in 
which the nodes are stacked in vertical columns of four 
(Fig. 20), is still used to obtain the solution, but is 
modified for plotting purposes. The z-coordinates are 
effectively restored by transformation to those describing a 
flat profile and the system is then rotated through the 
Slope angle. The columns of nodes are, subsequently, no 
longer vertical (Fig. 24). This change was instituted so 
that the entire plotting space on the AED 512 colour 
graphics terminal could be utilized. If the original 
coordinate system is used for plotting, then the available 
plotting area is much reduced. In all plots, the pronounced 
vertical exaggeration results in some distortion of 
isolines, particularly near lateral boundaries. 

The isopotential lines for the sloping system (Fig. 
24a) indicate that flow patterns are quite different from 
those in the horizontal system (Fig. 23a). Although water 
Movement is directed primarily toward the surface in each 
case, flow in the sloping profile has a downward component 
as well. Throughout the sloping profile water moves in 
response to gravitational potential; near the surface and 


downslope face, it moves also in response to forces produced 
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FORESTED SLOPE SIMULATION 
Evaporation : 2-5mm/day 
Transpiration : 2:5mm/day 


1000 2000 | 3000 
horizontal distance (cm) 
(a) TOTAL POTENTIAL (cm) 


Vertical exag: x9 


b 
a ° 


1000 2000 3000 
horizontal distance (cm) 


b) VOLUMETRIC WATER CONTENT (cm/cm? 


Figure 24. Evapotranspiration from a sloping profile -— total 
potential (a) and water content (b) -after 240 hr. 
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by evaporative demand. 

The sharp breaks in the isopotential lines at a depth 
of about 110 cm below the surface marks the limit of the 
surface evaporation and transpiration effects. Below this 
depth, flow toward the downslope face is indicated. At the 
bottom of the profile (0 to 30 cm) flow is toward the basal 
boundary. 

Soil water depletion patterns for the sloping (Fig. 
24b) and horizontal systems (Fig. 23b) are very similar. In 
both situations, the isolines are evenly spaced and parallel 
to the surface, indicating that water content increases 
uniformly with depth. At about 110 cm below the surface, it 
remains unchanged at 26 percent. The only real differences 
in water content between the horizontal and inclined systems 
result from evaporation at the downslope face of the sloping 
profile. Water content was reduced by about 10 percent at 
the node where the downslope face joins the surface 
boundary. 

Evapotranspiration from the same profile over a 20-day 
period, under the same conditions, was also simulated. 
Results indicate (Fig. 25) that the effects of evaporation 
after 20 days are still confined to the upper 110 cm of 
soil. However, movement of water above this depth is ina 
direction almost perpendicular toward the surface (eagr 
25a). Water content has been reduced by an additional 7 
percent at the surface, and by about another 5 percent at 


depth 60 cm below the surface (Figs. 24b and 255.) 
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FORESTED SLOPE SIMULATION _ 
Evaporation : 2-5mm/day 
Transpiration: 2-5mm/day 


1000 ~ 2000 
horizontal distance (em) 


(a) TOTAL POTENTIAL (cm) 


Vertical exag : x9 


A 


1000 2000 3000 
horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (cm/cm’ ) 


Figure 25. Evapotranspiration from a sloping profile - total 
potential (a) and water content (b) after 480 hr. 
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In the next four sets of experiments, the effects of 
locating trees on portions of the slope only will be 
examined. The soil, initial conditions, boundary conditions, 
time step size, and simulation period will all remain the 
Same as before. An evaporation rate of 2.5 mm/day will again 
be assumed for the entire slope. The only changing input 
variables will be the number and placement of activated sink 
nodes. 

Trees on the lower slope are simulated by activating 
only those sink nodes located below the midslope position. 
The 10-day simulation in this instance produced quite 
distinctive features in both the upslope and downslope 
Begions Of the profile (Fig. (26). Upslope, the effects’ of 
evaporation are evident down to a depth 70 cm below the 
surface. Downslope, the influence of evapotranspiration is 
manifest down to 110 cm below the surface. 

The equipotential lines signify that, uphill from the 
midslope position, water is moving upward in the top 70 cm 
of soil. Below this level, it is moving downslope (Fig. 
26a). On the downhill side of the midslope position, water 
is moving upward in the top 110 cm of soil, and moving 
downslope below. 

After 10 days, the reduction in soil water content is 
greater downslope than it is upslope (Fig. 26b). At a depth 
of 60 cm below the surface, for example, the upslope water 
content remains virtually unchanged from the initial 


condition of 26 percent. The downslope water content at this 
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TREES ON LOWER SLOPE 
Evaporation: 2.5mm /day 
Transpiration : 2-5mm/day 


1000, _ 2000 3000 
horizontal distance (cm) 


(a) TOTAL POTENTIAL (cm) 


Vertical exag: x9 


1000 2000 3000 
horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (cm/cm’) 


Figure 26. Trees on lower slope - total potential (a) and water 
content (b) after 240 hr. 
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depth, however, has been reduced by 4 percent to 22 percent. 

When data for this Pmeration (Fig. 26) are compared 
with those for the completely forested slope (Fig. 24), the 
most noticeable difference between them is the step-like 
change which occurs at the midslope position when trees are 
on the lower slope only. 

Results from the corresponding 20-day simulation (Fig. 
27a) show that, after 20 days, flow in the upper the portion 
of the profile is almost perpendicular to the slope. The 
regions affected by downslope evapotranspiration and upslope 
evaporation (Fig. 27) remain essentially the same as for the 
tO-day simulation (Fig 26).. 

Over the second 10-day period, volumetric water content 
at the surface has been reduced by an additional 4 to 6 
percent (Fig. -27b). The reduction downslope is greater than 
upslope. Above the midslope position, at a depth of 60 cm 
below the surface, reduction in water content is less than 
one percent. Downslope, at thd Same depth, the decrease is 4 
percent. 

For the next simulation, trees are located on the upper 
slope only. This effect is achieved by activating the Sink 
nodes above the midslope position, and de-activating them 
elsewhere. It is evident from the total potential field 
(Fig. 28a) for this system, that the response is analogous 
to that of the previous simulation, in which trees were 


located on the lower slope. 
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TREES ON LOWER SLOPE 
Evaporation: 2-5mm/day 
Transpiration : 2-5mm day 


x 


1000 es 2000 3000 
horizontal distance (cm) 


(a) TOTAL POTENTIAL (cm) 


Vertical exag: x9 


1000 “2000 3000 
horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (cm°/cm*) 


Trees on lower slope - total potential (a) and water 
content (b) after 480 hr. 
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TREES ON UPPER SLOPE 
Evaporation: 2:5mm/day 
Transpiration: 25mm /day 


1000. _ 2000 
horizontal distance (cm) 


(a)TOTAL POTENTIAL (cm) 


Vertical exag: x9 


ipod 1 dj gern Ah) 3000 
i t istance (c 
(b) VOLUMETRIC WATER CONTENT (cm'/cm?* ) 


Figure 28. Trees on upper slope - total potential (a) and water 
content (b) after 240 hr. 
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The effects of evapotranspiration extend to about 110 
cm below the surface, while downslope evaporation affects 
only the top 70 cm of soil. In these regions flow is 
directed upward toward the surface. Below these affected 
regions, flow is directed downslope (Fig. 28a). 

Soil water content, after 10 days, is less uphill from 
the midslope position than it is downhill. The difference is 
about 2 or 3 percent (Fig. 28b). At the surface, water 
content has been reduced by about 7.5 percent at upslope 
locations, and by about 5.5 percent downslope. Lowest water 
content (0.165) occurs at the point where the surface 
boundary joins the lower lateral boundary. 

The step-like change in data at the midslope position 
which was observed in the previous simulation (Figs. 26 and 
27) is also present in this simulation (Fig. 28). In this 
instance, however, the step is reversed. 

The isopotential lines for the 20-day simulation (Fig. 
29a) indicate that the soil regions affected by evaporation 
and evapotranspiration are the same as for the 10-day 
Simulation (Fig. 28a). They indicate that flow in these 
regions is upward and perpendicular to the slope; elsewhere 
flow is directed downslope. 

During the second 10-day period, soil water content at 
the surface is reduced by 6.5 percent in the 
evapotranspiration zone and by 4.5 percent in the 
evaporation zone (Fig. 29b). The net effect is a 4 percent 


difference in surface moisture content, between the two 
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TREES ON UPPER SLOPE 
Evaporation: 2-5mm/day 
Transpiration: 2:5mm/day 


1000 2000 
horizontal distance (cm) 


(a) TOTAL POTENTIAL (cm) 
Vertical exag: x9 


D9 
SS 
oo" 


1000 2000 3000 
horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (cm/cn?) 


Figure 29. Trees on upper slope - total potential (a) and water 
content. (b) atter 480 hr. 
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zones. At depth 60 cm below the surface the difference is 7 
percent. 

When all the sinks for the slope are de-activated, a 
"clearcut" condition is created and transpiration no longer 
occurs. The results from this simulation, shown in Fig. 30, 
are Similar in some respects to those obtained for the flat 
profile evaporation simulation (Fig. 21). In each case, 
after, 10 days, the effects of evaporation have reached a 
depth of 70 cm below the surface; above this level flow is 
directed upward toward the surface. For the clearcut 
condition on the sloping profile, flow below the 70 cm depth 
is directed downslope. 

After 10 days, soil water content at the surface has 
decreased by about 5 percent (Fig. 30b). For most of the 
profile, it increases uniformly with depth to about 70 cm 
below the surface. Below this depth, the water content 
remains unchanged. 

Results for the 20-day simulation (Fig. 31) show that, 
after 20 days, the effects of evaporation have not extended 
below the 70 cm depth. Above this depth, flow is 
perpendicular (upward) to the slope. Below, flow is still 
directed downslope (Fig. 31a). The surface soil water 
content has been reduced further by about 4 percent, 
bringing the total reduction from initial conditions at the 
surface to 10 percent (Fig. 31b). A decline in soil water 
content below the surface soil layer is also apparent, but 


it decreases with depth. 
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CLEARCUT SIMULATION 
Evaporation: 2:5mm/day 


Aa 


1000 1 1e2000 
horizontal distance (cm) 
(a)TOTAL POTENTIAL (cm) 


Vertical exag: x9 

Ay 
© oe. 
9 


3000 


1000 2000 
horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (cmi/cm’ ) 


Clearcut simulation - total potential (a) and water 


Figure 30. 
content (b) after 240 hr. 
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CLEARCUT SIMULATION 
Evaporation: 2:-5mm/day 


an aabrnas < 
2000 3000 


1000 
horizontal distance (cm) 
(a)TOTAL POTENTIAL (cm) 
Vertical exag: x9 


A 
Cs” 


cA 


1000 2000 3000 
horizontal distance (cm) 


“(b)VOLUMETRIC WATER CONTENT (cm’/om?) 


Figure 31. Clearcut simulation - total potential (a) and. water 
content (b) after 480 hr. 
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The preceding four sets of simulations serve to 
illustrate the possible effects of three different logging 
patterns on soil water flow and storage in a forested 
hillslope. The undisturbed forest (Figs. 24 and 25) can be 
compared with logging on the upper slope (Figs. 26 and 27), 
logging on the lower slope (Figs. 28 and 29), and with 
clearing of the entire slope (Figs. 30 and 31). The 
Simulations are very simplistic in that they contain fixed 
evaporation and transpiration rates for relatively long 
periods of time. No attempt is made, either, to accommodate 
the increase in ground surface evaporation which usually 
occurs when an area is cleared of trees. 

In spite of these limitations, the response for each 
Simulation is distinct, and reflects in a reasonable way, 
what happens in the corresponding field situation. If soil 
water content for clearcut conditions (Fig. 30b) is compared 
with soil water content for the completely forested slope 
Meigs 124b),,.it is clear that) removing:the trees has,resulted 
in a reduction in the drain on soil water. Under forested 
conditions, soil water content is diminished down to a depth 
of 110 cm below the surface. Under clearcut conditions, only 
the top 70 cm of soil is affected. In this region, at 
corresponding depths in the profile, soil water content:is 
higher for the clearcut than for the forested slope. These 
Simulations confirm, what has often been observed in the 
field, that removing trees results in an increase in soil 


water content. In some instances, especially where high 
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water-tables are present, the increase may be so great that 
it is difficult to restock the cleared sites with new trees. 

Patch logging and strip-cutting are two methods which 
have been used frequently in the past to provide a suitable 
environment for growing the next generation of trees at the 
Same time as mature trees are removed. Such methods are 
characterized by patches or strips of cleared or open areas 
alternating with strips or patches of trees. 

In order to simulate effectively the influence of these 
logging methods on soil water, the slope length used 4a 
previous simulations must be increased. By increasing the 
length to 100 m, two cut strips, each 20 m wide running 
across the slope and alternating with three forested strips 
also 20 m wide, can be accommodated. This configuration 
could be adapted quite easily to simulate the 
"wall-and-step" forest illustrated in Fig. 8. Except for the 
mesh size, all other conditions remain the same as before. 
The finite element mesh for this simulation contains 204 
nodes and 300 elements. 

The strip locations and the effects of this 
configuration on flow patterns are displayed in Pig. 2ar 
Below 120 cm from the surface, flow is everywhere parallel 
to the slope. -In the top 60 cm, flow is directed toward the 
surface. In this region, under the open strips, flow is 
almost perpendicular (upward) to the slope, whereas under 
forest it is inclined at an angle to the surface. Between 


depths 60 and 120 cm, two fairly distinct flow patterns are 
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PATCHCUT SIMULATION 
Evaporation : 2-5mm /day 
Transpiration : 2.5mm /day 


Rereontal aetonee (m) : 
(a ) TOTAL POTENTIAL (cm 


Vertical exag:x30 


20 40 , 60 80 100 
horizontal distance (m tee 
ENT (cm/cm’ ) 


Oo 
(b) VOLUMETRIC WATER CONT 


Figure 32. Patchcut simulation (two cut strips and three forested 
strips - total potential (a) and water content (b) 
after 240 hr. 
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evident. Under the open strips, water runs parallel to the 
Slope, but tends to move upward beneath the trees. This 
result 1S consistent with the fact that trees have a greater 
effect on water movement at depth, than does evaporation. 

Soil water depletion patterns for the patchcut 
configuration are quite pronounced (Fig. 32b), and are 
consistent with results obtained for previous simulations. 
Under the trees, soil water is depleted, to some extent down 
to a depth of 120 cm below the surface. For the cut strips, 
only the top 60 cm is affected (by evaporation). At the 
Surface, soil water content is 2 percent greater for the 
Open than for the forested strips. The reduction in water 
content decreases with depth for both sites. At depth 60 cm 
below the surface, soil water content 1S 26 percent under 
the open and about 23 percent under the forested strips. 

Comparison of the plots shown in Figure 32 with the 
Original data show that, although they are accurate over 
most of the region, the plots are not accurate in the 
vicinity of the lateral boundaries. The discrepency can be 
attributed to plot distortion caused by the large vertical 
exaggeration. 

If the positions of the open and forested strips are 
exchanged so that there are two forested and three open 
Strips, then the flow patterns and water content fields are 
also exchanged (Figs. 32 and 33). The flow patterns under 
trees in Fig. 32a are the same as those under trees in Fig. 


33a, even though the trees are on different parts of the 
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| PATCHCUT SIMULATION 
Evaporation : 2:'5mm/day 
Transpiration: 2:5mm/day 


40 60 80 
horizontal distance (m) 


(a)TOTAL POTENTIAL (c m) 


0) 100 


40 , 60 80 
horizontal distance(m) 


(b) VOLUMETRIC WATER CONTENT (cm/cm?) 


Figure 33. Patchcut simulation (three cut strips and two forested 
strips) - total potential (a) and water content (b) after 
ZeOentes 
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Slope. The same holds true for the open Strips. The 
Similarity of effects is particularly striking when water 
content fields are compared (Figs. 32b and 33b). The 
isolines, indicating points of equal soil water content, 
show identical patterns under the trees, where they are 
concave up, and under the open strips, where they are 
concave down. 

The flow systems examined 50 far are essentially 
isotropic, homogeneous, unsaturated systems. During these 
Simulations, problems relating to convergence or numerical 
Stability were not encountered. Usually, less than 10 
iterations were necessary to obtain convergence. The 
greatest number of iterations were realized during the first 
few time steps; for later time steps, fewer iterations were 


required. 
B. Combined Saturated-Unsaturated Flow 


Rainfall simulation 

The next simulation involves seepage resulting from 
rainfall, and introduces the effects of saturated flow into 
the formerly unsaturated flow system. The Yolo Light Clay 
Soil and the 64-node, 90-element mesh (Fig. 20) of earlier 
Slope simulations are used. No-flow boundaries are the same 
as before, and a positive flux equivalent to a rainfall rate 
of one cm/day is applied to each boundary node, with the 


exception of the corner nodes, along the upper surface. This 
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rate is slightly less than the saturated hydraulic 
conductivity, or drainage rate of the soil - a situation 
which implies that overland flow will not occur. The initial 
pressure potential for this problem is set everywhere equal 
to -100 cm (@=0.428). A time step of 6 hr is used to 
Simulate a 10-day period of rainfall. 

TheSinitial simulation is for rainfall into a 
horizontal box of soil from which there is no outlet. This 
experiment allows us to determine the change in soil water 
storage, and to compare it with the total volume of water 
added. 

The simulation output (Fig. 34a) indicates that, after 
10 days, total potential is highest at the surface, and 
lowest at the bottom, of the profile. Flow throughout most 
of the section is vertically downward. Near the lateral 
boundaries, flow is directed toward these boundaries as 
well. 

After 10 days, 10 em of rain has fallen onto the 
Mectile.eTheventire upper 120 cm of soil, with the exception 
of the regions near the lateral boundaries, is saturated 
(6=0.496) (Fig. 34b). This represents an increase in soil 
water content of about 7 percent, or an addition of 8.4 cm 
(0.07x120 cm) of water to the top 120 cm of soil over 10 
days. The balance, or 1.6 cm of added water, is contained in 
the bottom 60 cm of soil. These results are consistent with 


conditions that one might encounter in the field, following 


a l0-day rainfall. 
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CLEARCUT SIMULATION 
Rainfall : ]Omm/day 


HEIGHT (cm) 


1000 2000 3000 
horizontal distance (cm) 


TOTAL POTENTIAL (cm) 


Vertical exag : x8 


HEIGHT (cm) 


O 
horizontal distance (cm) 


VOLUMETRIC WATER CONTENT (cm*/cm?) 


Rainfall on a flat profile - total potential (a) and 


Figure 34, 
water content (b) after 240 hr. 
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The 64-node, 90-element mesh for a Sloping section is 
employed to simulate outflow from the seepage face. In this 
simple problem, nodes numbered two and three, located along 
the downslope lateral boundary, are seepage face nodes. The 
only other difference between this simulation and the 
preceding one is that, in this problem, rainfall is applied 
to the corner node (#4) at the intersection between the 
lower lateral boundary and the upper surface, in addition to 
the other surface nodes. 

The total potential field that exists at the end of the 
10-day period (Fig. 35a) indicates that flow is directed 
downslope toward the lower lateral boundary. The top 120 cm 
of soil is saturated (6=0.496) except near the lateral 
boundaries (Fig. 35b). It is not clear why there 1s a pocket 
of unsaturated soil below seepage face node #3 at the 60 cm 
depth. One would expect the "Saturated" isoline to continue 
along its path parallel with the slope, and terminate at a 
depth of about 120 cm on the lower lateral boundary. This 
assumption is supported by evidence from Fig. 35a, where the 
isopotential lines show that flow is directed toward the 
basal portion of the lower lateral boundary. 

Outflow did not commence until 4.5 days after the rain 
Started (Fig. 36a). The hydrograph shows an abrupt rise 
during the first two time steps after flow started, followed 
by a somewhat erratic, but generally steady, increase over 


the rest of the simulation period. 
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CLEARCUT SIMULATION 
Rainfall : 1Omm/day 


1000 2000 
horizontal distance (cm) 
(a) TOTAL POTENTIAL (cm) 


Vertical €xag: x9 


zx 
1000 2000 3000 


horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (cm?/cm?) 


Figure 35. Rainfall on a sloping profile - total potential (a) 
and water content (b) after 240 hr. 
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Since there is no physical explanation for the apparent 
anomaly of an unsaturated pocket below the saturated zone at 
the lower lateral boundary, additional simulations were run 
to investigate the problem further. It was hypothesized that 
the anomaly was related to the coarseness of the mesh and to 
insufficient seepage face nodes. Accordingly, the problem 
was rerun uSing a mesh consisting of 217 nodes and 360 
elements. In this case, the number of seepage face nodes was 
increased to five. 

The results show (Figs. 35a and 37a) that changing the 
mesh size has little effect on flow patterns - the total 
potential fields are almost identical. 

The change does affect the volumetric water content of 
the profile. Saturation extends to slightly below the 120 cm 
depth for the fine mesh (Fig. 37b) compared with saturation 
to slightly above the 120 cm depth for the coarse mesh (Fig. 
35b). Noticeable changes are apparent near the lateral 
boundaries. At the upslope boundary in the fine mesh 
simulation, saturation reaches a depth of about 100 cm (Fig. 
37b). In the coarse mesh simulation, saturation does not 
occur anywhere within 250 cm of the upper lateral boundary. 

At the lower lateral boundary, for the fine mesh, the 
unsaturated pocket is more prominent than before, and the 
seepage face extends to a depth of only 40 cm, compared with 
70 cm for the coarse mesh. 

When the hydrograph for the fine mesh Simulation (Fig. 


36b) is compared with that of the coarse mesh simulation 
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CLEARCUT SIMULATION 
Rainfall : 10mm day 


1000 2000 
horizontal distance (cm) 
(a)TOTAL POTENTIAL (cm) 


Vertical exag: x9 


ce 
3000 


OO 
horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (cm*/cm?) 


Figure 37. Rainfall on a sloping profile - total potential (a) and 
volumetric water content (b) after 240 hr. Obtained 
using a 217-node, 360-element mesh. 


ems 


ee 


a Se te —_ 


onan imentea niet 


oo. 
Crist aorarats 
(may. TANG 


My. eee p< a ee ae 


Our, 


=. sing Te 


a 


host 46) ipksaoeuact ns 


Serhasdt) > +) i ae 


2 


TALS can 


174 


(Fig. 36a), it is clear that the hydrograph is influenced by 
mesh size. Outflow from the seepage face begins two days 
earlier when the fine mesh is used. It does so because the 
first seepage face node is located at depth 30 cm, compared 
with 60 cm in the coarse mesh. The hydrograph for the fine 
mesh is also much smoother - it lacks the abrupt rise and 
irregularities that characterize the coarse mesh hydrograph. 

Maximum flow rate is greater for the coarse mesh 
Simulation (5.8 cm?/hr). It would appear from the foregoing 
results that more water is assigned to storage and less to 
seepage when a fine mesh is used. The fine mesh simulation 
did not provide any clues as to why an unsaturated region 
persists at the lower seepage face. 

Since the hydraulic conductivity of the soil is 
governed by soil wetness, the volume of outflow and the time 
at which seepage commences must obviously be controlled by 
the initial condition of the soil, as well as by the amount 
@e inflow. ‘This is particularly true in the vicinity of the 
seepage face. The effects of changing the initial conditions 
at the seepage face were examined in an ensuing fine mesh 
simulation. Only one change was implemented - the initial 
pressure potentials at the seepage face, formerly set at 
-100 cm, were reinitialized to -50 cm. 

This change in initial conditions produced only a minor 
change in the response of the system, a change which was 
confined to the region near the seepage face. It did not 


affect the unsaturated pocket. Outflow began a day earlier 
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(Fig. 36c) than before. The hydrograph is slightly higher 
than, and parallels that of, the previous simulation (Figs. 
36b and 36c). 

Another fine mesh simulation, in which wetter initial 
conditions at the seepage face were assumed, was run. 
Initial pressure potentials were set to -10 cm at the 
seepage face, and to -50 cm at the adjacent nodes. 

After 10 days, the modification produced only a slight 
change in the total potential field (Figs. 37a and 38a), but 
it was sufficient to cause complete saturation at the 
seepage face (Fig. 38b) and increased outflow (Fig. 36d). 
Only the total potential (Figs. 37a and 38a) and water 
content (Figs. 37b and 38b) in the vicinity of the seepage 
face were changed. 

Initially, the hydrograph (Fig. 36d) is similiar to 
those for the other fine mesh simulations (Figs. 36b and 
36c). Outflow begins at the one half day mark and, for 34 
time steps, takes place through the uppermost portion of the 
Seepage face. After that time, the entire lower lateral 
boundary rapidly becomes saturated, resulting in a 
pronounced increase in outflow from the entire seepage face 
(Fig. 36d). When these wetter conditions prevail, more 
iterations are required to obtain convergence. 

The simulation was extended beyond the 10-day period, 
but it failed at time step 42 (day 11). Failure was 


characterized by oscillations which led ultimately to 


divergence. 
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CLEARCUT SIMULATION 
Rainfall: IOmm/day 


xz 


1000 2000 
horizontal distance (cm) 


(a) TOTAL POTENTIAL (cm) 


Vertical exag:x9 


1000 2000 3000 
horizontal distance (cm) 


(b) VOLUMETRIC WATER CONTENT (€m*cm?) 


Bigure 38. Rainfall on a sloping profile with pj = -10, -50, and 
-100 cm - total potential (a) and volumetric water 
content (b) after 240 hr. Fine mesh simulation. 
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Changing boundary conditions 

Results from the rainfall simulations suggest that 
SUBFEM is sensitive to mesh size and to initial conditions. 
Results from additional, exploratory test simulations 
Suggest that the model is also highly sensitive to changing 
boundary conditions. Until now, the boundary conditions for 
each simulation remained unchanged for the entire simulation 
period. This constraint allowed the behaviour of each 
hydrologic process to be examined in isolation. In the 
following simulations the effects of changing boundary 
conditions are investigated. 

The first run involving changing boundary conditions 
Simulated 3 days of rainfall followed by 10 days of 
evaporation. Initial conditions corresponded to those that 
existed after 7 days, when the seepage face approached 
Saturation, in the previous simulation. The simulation 
proceeded smoothly until the new boundary conditions were 
introduced, at which time the simulation failed. Oscillation 
occurred and no solution was obtained. 

At first, it was assumed that the problem was related 
to time step size. However, when the time step size was 
adjusted from 6 hr to 3, 0.5,and 0.1 hr, there was no 
improvement. It was then hypothesized that either there was 
an error in the program or that the changes in boundary 
conditions were too abrupt to allow a solution to be 
obtained. The program was carefully checked for errors but 


none were detected. Consequently, it was decided that a 
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detailed investigation of changing boundary conditions was 
warranted. 

Five simulations were run for each of two changing 
boundary conditions: the first in which rainfall changed 
trom 10 mm to 2 mm/day, and the second in which rainfall 
changed from 10 mm to 0 mm/day. Initial conditions were set 
at the conditions which existed at the end of the 10-day 
rainfall simulation for the Yolo Light Clay and fine mesh 
(Figs. 36d and 38). Two time periods were examined - the 
first hour, and the first day following the time at which 
the boundary conditions changed. Time step sizes of 6, 3, 1, 
0.1, and 0.01 hr were used. 

When the rainfall rate was changed from 10 mm to 2 
mm/day, there was a marked drop in outflow (Fig. 39) during 
the first time step (At=0.01 hr). There was no appreciable 
change in the water content field, but there was a 
Significant change in pressure potentials in the saturated 
zone, particularly at the surface nodes. During the first 
time step after the boundary conditions were changed, 
positive pressure potentials at the surface nodes were 
reduced from about 120 cm to 35 cm. This caused the isolines 
for total potential to intersect the upper surface at an 
angle steeper than before. 

The effects of time step size on the Simulation are 
shown in Figs. 39 and 40. It is clear that, during the 
latter part of the simulation, the same solution is obtained 


regardless of time step size. However, during the early part 
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of the simulation, immediately following the change in 
boundary conditions, time step size is important. An 
improved representation of the recession curve is obtained 
each time the step size is reduced. 

The results obtained for the simulations in which 
rainfall ceased after 10 days (Figs. 41 and 42) are similar 
to those for simulations in which rainfall was reduced from 
10 mm to 2 mm/day after 10 days (Figs. 39 and 40). Each 
hydrograph shows a pronounced drop immediately following the 
change in boundary conditions. Where rainfall persists, 
outflow continues beyond the first day (Fig. 40). When the 
rain stops, outflow ceases 13 to 18 hr later (Fig. 42). 

At the end of the first time step following cessation 
of rainfall, the pressure potential in much of the saturated 
region is uniform, with a value of 14 cm. This is 
Significantly lower than most of the positive pressure 
potentials which existed at the beginning of the time step. 
The drop in pressure potential causes the total isopotential 
lines in the upper 120 cm of soil to bend toward the 
perpendicular relative to the surface boundary (Figs. 38a 
and 43). The water content field is not significantly 
affected during the first time step. 

Time step size-has the same effect on this simulation 
as it had on the simulation in which rainfall was reduced 
from 10 mm to 2 mm/day. Results differ for time steps 
immediately following changes in boundary conditions, but 


the differences diminish as the simulation proceeds (Figs. 
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Rainfall : 10mm to Omm 
At=0-01 


1000 2000 
horizontal distance (cm) 


TOTAL POTENTIAL (cm) 


Figure 43. Total potential at end of first time step 
following cessation of 10 cm rainfall. 
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41 and 42). Only two time steps were completed for the 
Simulation in which At=1 hr. The simulation was terminated 
after 50 iterations during the third time step, when 
convergence was very slow. 

Results from the foregoing simulations for combined 
Saturated-unsaturated flow illustrate the classical concept 
of interflow. As rainfall continues, a saturated layer 
develops over an unsaturated zone in the originally 
unsaturated soil. Eventually, outflow issues from the 
downslope lateral boundary. Although the situation lacks the 
conventional prerequisites, such as macropores (root 
channels and animal burrows), layered soil, and hardpan 
layer, to produce significant amounts of interflow, there is 
no question that interflow is occurring. The outflow comes 
directly from the upper saturated layer, and not from 
overland, nor via a groundwater system at greater depth. 

It is evident from the simulations completed so far for 
Saturated-unsaturated eee that SUBFEM is sensitive to mesh 
size, time step size, initial conditions, and changing 
boundary conditions. In the next section, the effects of 
changing some of these conditions, and the effects of 


changing saturated hydraulic conductivity will be examined 


in some detail. 
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Sensitivity analyses 


Time step size 

The problem selected for sensitivity analyses entails 
Simulation of 5 days of evapotranspiration (3 mm/day ) 
followed by 10 days of rainfall (10 mm/day). The 
evapotranspiration consists of 1.5 mm/day evaporation and 
1.5 mm/day transpiration. In the first set of simulations, 
the effects of different time steps sizes are investigated. 
For this purpose, the 64-node, 90-element mesh (Fig. 20) and 
Yolo Light Clay are used, and the initial pressure 
potentials are set everywhere equal to -100 cm (6@=0.428). 
Five time step sizes: 1, 3, 6, 12, and 24 hr, are simulated. 

When data for different time step sizes at the end of 
each day are compared, it is found that, in general, the 
pressure potential fields (and hence the total potential and 
water content fields) are very similar. Such differences 
that occur are usually of the order of about 2 or 3 cm, and 
appear after the first two or three days of simulation. 
Larger differences (10 to 15 cm), which occur at some nodes 
in the vicinity of the seepage face, ordinarily arise as the 
water content at the seepage face approaches saturation. 
These differences are sustained for some period following 
Saturation, but subsequently diminish with time. The 
comparison which showed the greatest differences was thaweor 
At=1 hr and At=24 hr. 

Discharge from the seepage face commenced 10 1/2 days 
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simulations of all time step sizes except, perhaps, for the 
At=24 hr simulation. It is clear that, with the exception of 
the At=24 hr curve, there is very little difference between 
plots for different time step sizes, beyond the beginning of 
day 11. However, it is equally evident that the rising limb 
of the hydrograph is strongly dependent on time step size - 


its slope increases as the step size decreases. 


Mesh size 

Four mesh configurations are used to determine the 
effects of mesh size on output. Two have already been 
referred to - the 64-node, 90-element mesh and the 217-node, 
360-element mesh. The first has a maximum of two, and the 
second, five seepage face nodes. Two other, intermediate, 
mesh patterns are introduced; they consist of 105 nodes on 
160 elements, and 114 nodes on 172 elements. These two 
meshes are identical over much of the region being 
Simulated. They differ only in the vicinity of the seepage 
face, where the 105-node mesh has 3 seepage face nodes. In 
contrast, the 114-node mesh has 7 seepage face nodes. The 
simulations are done for At=3 hr and At=24 hr. 

The results from this set of simulations indicate that, 
after 15 days, for each mesh configuration, the total 
potential field obtained when At=3 hr is almost identical to 
that obtained when At=24 hr. Differences at any given 


location in the simulated region are rarely greater than 3 
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When results for the same time step size and different 
mesh sizes are compared, some differences in potential 
fields are apparent, most noticeably in the proximity of the 
seepage face. Differences in potential tend to be greater at 
the bottom of the region than at the top. 

A point which is common to all mesh configurations 
tends to have similar total potential values for different 
meshes. The isopotential lines are different, however, 
because each mesh contains a different number of nodes or 
data points. As this number increases, a greater degree of 
interpolation within the region is possible. The result is a 
more accurate representation of the isopotential lines. 
Consequently, the mesh with the greatest nodal density (the 
217-node mesh) should provide the most accurate solution for 
a given time ’step size. 

The two intermediate-sized meshes give identical 
results for the region 600 cm from the seepage face and 
beyond. In the region between the seepage face and the 600 
cm horizontal distance mark, differences in total potential 
between the two meshes tend to increase as the distance from 
the seepage face diminishes. 

The time at which outflow commences from the seepage 
face is governed by the proximity to the surface of the 


first seepage face node. The depths to the first seepage 


face node for the 64-, 105-, 114-, and 217-node systems are 


60, 45, 22.5, and 30 cm respectively. 
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Results from the simulations show (Figs. 45 and 46) 
that outflow begins soonest when the 114-node simulation is 
run, and latest when the 64-node system is used. Outflow 
from each mesh commences in the order in which the first 
seepage face node is tapped. The hydrographs for At=3 hr are 
qualitatively similar to those for At=24 hr. 

The data show (Figs. 45 and 46) that the earlier 
outflow starts, the more gradual its increase afterwards. In 
contrast, the longer the delay before discharge begins, the 
more abrupt is the initial increase in outflow likely to be. 

These observations can be related to the fact that if 
water cannot be discharged from the system, then it must be 
Bored. Thus, in the 64-node system, water will not be 
discharged until the first seepage face node (at depth 60 
cm) reaches saturation. When this condition is attained, 
some of the water stored in the saturated layer above is 
released immediately. When the 114-node system is used, the 
Saturated layer has to be only 22.5 cm thick before outflow 
begins. In this case, less water is discharged than when the 


Saturated layer is 60 cm thick. 


Initial conditions 

The 114-node, 172-element mesh with 7 seepage face 
nodes was chosen to test model sensitivity to changes in 
initial conditions. In this analysis, a 3-hr time step is 
used and Pane initial conditions are compared. At the 


beginning of the first simulation, the pressure potential is 


set everywhere equal to -100 cm (0=0.428). 
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In the other simulations, the initial pressure 
potential fields are not uniform. Instead, each of them 
corresponds to an upper moist region 45 cm deep overlying a 
dryer region 135 cm deep. The moist region extends downward 
along the seepage face, so that there is a uniformly moist 
band extending 75 cm behind the face. Within both moist and 
dryer regions initial pressure potential is uniform. For 
three separate simulations, the initial pressure potential 
is set at -150 cm in the dryer region, and at -5, -10, and 
-50 cm in the moist region. 

The results from the simulations indicate (Fig. 47) 
that both the time at which outflow commences and the volume 
of flow are affected when initial conditions are changed. 
For a given set of conditions, it is evident that the 
greater the degree of surface saturation imposed at the 
beginning of the simulation, the earlier seepage begins. 
Thus, when the initial pressure potential is set at -5.0 cm 
in the upper region and -150 cm in the lower region (yWj,=-5 
cm/-150 cm), the discharge from the seepage face begins 
during day 5 of the simulation. When the initial pressure 
potential is prescribed uniformly as -100 cm, outflow does 
not begin until day 7. It also appears that greater initial 
flows, and greater flows generally, are linked to early 
initiation of outflow. There is practically no hai F AIS 
between hydrographs obtained when w;=-10 cm/-150 cm and 


wW,=-5cem/150 cm [c) and d) in Fig. 47]. 
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Hydraulic conductivity 

The effects of changing hydraulic conductivity were 
evaluated using a 3~hr time step, and the 114-node, 
172-element mesh with 7 seepage face nodes. Initial pressure 
potentials were set at -50 cm along the seepage face and 
upper two rows of nodes, and at -150 cm elsewhere. 

Eight simulations were completed. The first concerns 
Yolo Light Clay which has a saturated hydraulic conductivity 
of 0.04176 cm/hr (slow). This is the same simulation which 
was used in the previous experiment for initial conditions 
analysis (Fig. 47b). It will now serve as the reference 
Simulation for comparing the effects of changes in 
conductivity. The hydrograph is repeated (Fig. 48a) for this 
purpose. 

In the second simulation, the saturated hydraulic 
conductivity is reduced by one order of magnitude to 
0.004176 cm/hr (very slow). For this problem, the drainage 
rate is now less than the rainfall rate - a condition which 
gives rise to overland flow. 

The reduction in conductivity causes less water to 
reach the lower regions of the profile (Fig. 49b) than 
previously (Fig. 49a), and very little penetrates beyond 
depth 135 cm. In both simulations, a saturated (6=0.496) 
layer together with positive pressure potentials develop in 
the upper portion of the profile. However, the positive 
pressure potentials for the simulation with the lower 


saturated conductivity are much higher. This is a measure of 
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Figure 49. Water content in Yolo Light Clay after 10 days 
of rainfall, for two values of saturated 
hydraulic conductivity (K,). 
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the model's response when no provision is made for overland 
flow. Since rainfall in excess of the drainage rate is not 
removed through an overland flow routine, the model tries to 
accommodate the excess water in some other way - evidently 
by generating high positive pressure potentials. 

ane=results®also*indicate that /"for “atgiven Sseerot 
conditions, a reduction in saturated conductivity leads to a 
thinner saturated layer (Figs. 49a and b). Because of this, 
and since less water is moving downward through the lower 
regions of the profile, one would expect more water to 
appear at the seepage face. Inspection of the hydrographs 
for the two simulations (Figs. 48a and b) shows that this 
is, indeed, the case. 

For the third simulation, the reference conductivity is 
doubled to 0.08352 cm/hr (slow). When the results are 
compared with output from the reference simulation, two main 
trends emerge. First, more water is stored in the profile, 
especially in the basal portion, and second, outflow from it 
through the seepage face is reduced (Figs. 48a and c). 
Because the transmitting capability of the soil has 
increased, water has more rapid access to the unfilled 
storage capacity in the lower portion of the profile. 
Consequently, less water issues from the seepage face. The 
Saturated layer extends from the surface to a depth of about 
7OFCmM. 

When the saturated conductivity is increased one order 


of magnitude above the reference conductivity to 0.4176 
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em/hr (moderate), the entire amount of rainfall is used to 
satisfy pea aa b ie Soil water storage. In this case, 
Saturation is not reached anywhere and outflow does not 
eecur . 

Similar results are obtained when a simulation is 
performed on data pertaining to reconstituted Halewood Sandy 
Loam (Scholl and Hibbert, 1973). The saturated conductivity 
for this material is about 0.45 cm/hr (moderate). 

Some outflow is obtained when the problem is run using 
a Saturated conductivity of 0.20 cm/hr (slow). However, only 
a small region, in the vicinity of the seepage face, reaches 
Saturation. Consequently, outflow is minimal and does not 
begin until near the end of the simulation (Fig. 48d). 

The remaining two conductivity Simulations show that, 
if inflow is adequate, large saturated conductivity values 
-result in steady, saturated flow. The first of these 
Simulations concerns Calvin Silt Loam (Troendle, 1970) which 
has a saturated conductivity of 20 cm/hr (rapid). The soil 
used in the other simulation is Bodine Silt Loam (Huff et 
al, 1977) which has a saturated conductivity of 60 cm/hr 
(very rapid). 

In both simulations, saturation occurs before they are 
completed. For the Calvin soil, saturation is attained 
during day 11 (Fig. 48e), and for the Bodine soil, during 
day 13 (Fig. 48£). When saturation is reached, all pressure 
potentials are zero or positive, and total outflow equals 


total inflow. This condition is maintained for the remainder 
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of the simulation (Figs. 48e and f), 

The time at which outflow from each soil begins is 
related to the available storage capacity of the soil. This 
capacity, in turn, depends on the soil's initial water 
content and its total porosity. Data used in the Bodine soil 
Simulation indicate that, at the start of the simulation, 
this soil has the capacity to absorb about 6.5 cm of 
additional water. The corresponding value for the Calvin 
soil is about 9.5 cm. Because water is being applied and 
absorbed at a rate of 1 cm/day from day 6 on, one would 
expect outflow to occur (at the inflow rate) from the Calvin 
soil 6.5 days later, and from the Bodine soil, 9.5 days 
later. The simulation results (Figs. 48e and f) are 


consistent with these expectations. 


C. Model Versus Prototype 

Attempts were made to simulate a field experiment which 
waS carried out on the Fernow Experimental Forest, West 
Virginia, during the late sixties. The field experiment 
(prototype) lends itself to simulation using SUBFEM because 
it is two-dimensional and because data from it are quite 
comprehensive. 

The field study was conducted on a fairly uniform 
north-facing slope with a gradient of about 19°. The soil is 
predominantly deep Calvin Silt Loam, with inclusions of 
Dekalb loam, overlying Catskill sandstone and shale. Forest 


vegetation consists of oaks, sugar maple, yellow poplar, 
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black cherry, and beech. The plot extends from the ridge to 
a horizontal distance downslope of 53.3 m. At the lower end 
of the plot, a collection trench 15.2 m wide and about 1.8 m 
deep was excavated along the contour perpendicular to the 
Slope, to intercept soil moisture draining downslope 
(Troendle, 1970). 

Backfilling and sealing of the trench was designed so 
that a wall of pea-sized gravel impinged on the upslope face 
of the trench, and all water issuing from this face was 
collected in a trough at the bottom of the trench. The water 
was then transferred from the trough into an HS flume. 

To prevent lateral movement of water out of the plot, 
its lateral boundaries were trenched to bedrock and 
backfilled with a mixture of Bentonite clay and soil. 
Instrumentation was installed to measure soil water content 
and to monitor perched water tables. 

For the purpose of simulation, it was assumed that the 
basal boundary of the hillside coincided with the bedrock at 
depth 1.8 m. It was further assumed that this boundary ran 
parallel with the surface from the ridge to the trench. 
Another, vertical, impermeable boundary was presumed to 
exist at the topographic divide, or ridge. 

Soil descriptions for the plot indicate a multi-layered 
profile. For simulation purposes, however, only two layers 
are considered - the upper, which extends from the surface 
to a depth of 62.5 cm, and the lower, which lies between 


depths 62.5 and 175 cm. Since saturated conductivity values 
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for the upper horizons are not given, a value of 170 em/hr 
for the upper layer is assumed. This is an order of 
magnitude greater than the saturated conductivity for the 
lower layer. Each layer is considered to be isotropic and 
homogeneous. 

The simulated event concerned a rainstorm which 
Cecurred:’on 8-9 August,1969. During this storm 7.75 cm of 
rain fell, in less than 12 hr, onto an Originally dry soil. 
The simulation was accomplished using a 272-node, 
462-element mesh. 

The distribution of rainfall for the storm is given as; 
eo etc) 2 1.127 ,10R63', 90551; and >0.13 cm/hr. These rates, 
which are applied over 2-hr increments for a total of 12 hr, 
represent the boundary conditions of the problem. On August 
6, volumetric soil water content was about 18 percent in the 
lower layer and between 19 (upslope) and 24 (downslope) 
percent in the upper layer. The water contents define the 
initial conditions of the problem. A time period of 12 hr 
was Simulated using a time step of 0.1 hr. 

Of the total amount of precipitation that fell during 
the August storm, only 5.5 percent appeared as outflow from 
the plot, the rest recharged the soil (Troendle,1970). In 
the simulation , no point on the hillslope reached 
Saturation, and no outflow was obtained. 

A further simulation of the August storm was carried 
out. It differed from the preceding one in that initial 


water contents for the surface soil layer were prescribed as 
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26 and 28 percent, instead of 19 and 24 percent. 

The results from this simulation, also, indicated 
unsaturated conditions and a complete lack of outflow from 
the hillslope for the duration of the 12-hr period. 
Evidently, the entire 7.75 cm of rainfall was absorbed by 
the soil. Simulated water content data showed that, by the 
end of the 12-hr period, water content in both layers had 
increased by about 4.5 percent, or by a total of 
approximately 8.0 cm - a figure somewhat higher than the 
eriginal weoAS cm .of rain. 

The data for each time step show that, at the bottom of 
the profile, no appreciable change in water content occurred 
until after time step 60 (6 hr). During that time, at the 
Surface, the downslope pressure potential increased from 
-4944 to -371 cm, and the upslope pressure increased from 
-7210 to about -405 cm. Thus after 6 hr, soil water content 
at the surface was close to field capacity (wW=-339 cm). 

A fairly steady Ba Seenn was evident for the period 
between 6 and 10 hr. At depth 25 cm, pressure potential did 
not change very much - remaining at about -616 to -620 cm. 
Below depth 25 cm, soil water content increased at all 
depths. Above depth 25 cm, the pressure potential and water 
content decreased slightly. 

After 10 hr, when the rainfall rate was reduced to 0.13 
em/hr, the pressure potential at the surface decreased more 
rapidly, indicating that the upper soil layer was draining 


at a rate faster than the supply rate. 
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It is not clear why there was outflow from the 
prototype system, yet none from the model system. Some 
possible explanations are offered below: 

1. In the simulation, each soil layer was considered to be 
homogeneous, thus its spatial variability was not taken 

. into account. If the layers are not homogeneous, then 
the given point measurements of saturated hydraulic 
conductivity may not be applicable over the entire 
hillslope. 

2. Each soil layer was assumed to be isotropic for 
Simulation purposes. This may not be the case for the 
prototype. If the horizontal conductivity exceeds the 
vertical conductivity, water will tend to move downslope 
toward the seepage face faster than it will vertically. 

3. The model simulates only flow through porous media. It 
is possible that a different mechanism is operating to 
produce outflow from the prototype, e.g. turbulent flow 
through macropores such as root holes and animal 
burrows. 

The simulation just described was the largest attempted 
during this study. Running the problem using the 272-node, 
462-element mesh and 120 time steps entailed 8 minutes of 
CPU time and 1996 page-minutes of virtual memory. The input 


data files for this problem are displayed in Appendix D. 
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VI. SUMMARY AND CONCLUSIONS 


A. Summary 

A physically based, distributed, subsurface flow finite 
element model (SUBFEM), in which vegetation appears as a 
integral part, was developed to simulate the effects of 
different vegetation patterns on soil water distribution and 
Streamflow. Water withdrawal by trees was simulated by means 
of sinks located at appropriate near-surface nodes within 
the finite element mesh. Each sink may represent the 
activity of a single tree or groups of trees. In this study 
each sink represented water withdrawal by a group of trees. 

Data required to drive the model include section 
geometry and, for each soil or rock layer, saturated 
hydraulic conductivity and the relation between volumetric 
water content and pressure potential. Rainfall (or 
snowmelt), evaporation, and transpiration rates over time 
are also required. For each time step the program prints out 
the pressure potential field and boundary fluxes. Boundary 
fluxes occur at the upper boundary (ground surface) and at 
the seepage face. The program also prints out the pressure 
potential, the total potential, and the volumetric water 
content fields, together with the boundary fluxes, that 
exist at the end of the simulation. 

Initial tests on the model, for one-dimensional and 
two-dimensional infiltration into originally unsaturated 


soil, produced results that compared favourably with those 
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of other investigators who had simulated the same problems. 
These tests also successful lyademonst nated the concepts of 
advancing wetting front and infiltration Capacity. 

A series of simulations were completed to determine the 
effects of transpiration and evaporation, both separately 
and combined, on water distribution in a large box of soil. 
The simulations revealed that a smaller volume of soil is 
affected by evaporation than by transpiration for a fixed 
evaporation or transpiration demand. Potential gradients, on 
the other hand, are greater during evaporation. 

The next set of experiments concerned 
evapotranspiration, over 10- and 20-day periods, from 
Sloping profiles. In these simulations, trees were assumed 
to occupy different positions on the slope: lower slope 
only, upper slope only, completely forested, clearcut (no 
trees), and patchcut. 

Quite distinct patterns for both total potential and 
volumetric water content were obtained. They reflected the 
differential drain on soil water by trees and by 
evaporation. Soil water content was less, and the volume of 
soil affected greater, under trees than under open 
conditions. The isopotential lines indicated that, when 
evaporation or transpiration occurs, flow is directed up 


toward the surface in the upper portion of the profile, and 


downslope at lower depths. 
The first simulation in which saturated flow was 


considered involved rainfall into a horizontal box from 
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which there was no outlet. The increase in soil water 
content over the 10-day rainfall period closely approximated 
the total amount of rainfall. 

In a subsequent problem, rainfall over a 10-day period 
onto a sloping profile was simulated. This time outflow was 
permitted. Although the simulation involved porous media 
flow only, results from it illustrated the classical concept 
of interflow. A saturated layer developed over an 
unsaturated zone in an originally unsaturated soil, and 
seepage eventually emerged from the upper portion of the 
downslope face. . 

Simulation was used in an attempt to reproduce results 
from a field experiment conducted on the Fernow Experimental 
Forest, West Virginia, during the late sixties. A single 
Storm event,°in which 7.75 cm of rain fell onto an 
Originally dry soil in less than 12 hr, was Simulated. 

Field data showed that 5 percent of the precipitation 
appeared as outflow from the plot. The simulation data 
indicated that no outflow occurred during or following the 


storm, and that the entire contribution from the storm was 


retained in the soil. 


B. Conclusions 


The conclusions that follow are based on simulations of 


flow through isotropic, homogeneous porous media. Results 


from these simulations suggest that the primary usefulness 


of SUBFEM rests in its capability to simulate, qualitatively 
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at least, the processes of infiltration, evaporation, 
transpiration, and interflow. Because of this, and since it 
is two-dimensional, it can best be used to study hydrologic 
processes on vegetated or partially vegetated hillslopes. 
The simulation of the plot study on Fernow Experimental 
Forest exemplifies this type of application. 

The data requirements of the model are considered to be 
one of its main limitations. For example, transpiration 
rates for different tree species under different climatic 
conditions are difficult to obtain and, generally, are not 
available. Assigning suitable saturated conductivity values 
to soil strata is also a problem. 

Experience shows that SUBFEM can easily cope with 
large, unsaturated flow systems. For these simulations, few 
iterations per time step are required, and only rarely are 
problems encountered that relate to convergence or numerical 
instability. Some difficulties arise when 
Saturated-unsaturated flow systems are simulated - notably 
following certain changes in boundary conditions. 

Time step size, under certain conditions, is critical 
and will determine whether convergence to a solution is 
possible. 

Mesh coarseness can influence the time at which outflow 
from an originally unsaturated system begins. Sensitivity of 
response, in this case, is related to the number and 
position of nodes at the seepage face. For infiltration 


problems, the more nodes there are at the seepage face che 
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earlier outflow begins. 

The sensitivity of two model parameters were examined: 
initial conditions and saturated hydraulic conduc tivt tye tof 
the two, saturated hydraulic conductivity proved to be the 
most sensitive. 

Both the time at which outflow commences and the volume 
of flow from. an originally unsaturated system are affected 
by changes in initial conditions. If the region in the 
vicinity of the seepage face, and the upper portion of soil 
are wetter than the soil below, then together they serve as 
a conduit for additional water that enters the system. In 
general, the wetter the initial conditions, the earlier 
outflow begins. 

SUBFEM is very responsive to changes in saturated 
hydraulic conductivity. If the value of this parameter is 
high enough, and the porous material is homogeneous and 
isotropic, the entire profile may be recharged before 
outflow occurs - at either the saturated or inflow rate, 
whichever is lower. Afterwards, if the same inflow rate is 
maintained, steady state conditions prevail. 

Although the model is designed primarily to simulate 
forest site conditions under which soil absorbs rainfall at 


most intensities, it was established that, if the rainfall 


rate exceeds the saturated conductivity, provision must be 
made in the model to route excess water as overland flow. 


There are indications that provision should also be made to 


simulate flow processes other than porous media flow, 1.e. 
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turbulent flow through macropores such as root holes and 


animal burrows. 


C. Suggestions For Further Research 


ao 


SUBFEM should be subjected to more rigorous testing 
by comparing output for a given field situation with 
observed data from a well-instrumented site. 

The subsurface flow model may prove to be a useful 
tool for determining average saturated hydraulic 
conductivity values for hillslopes. 

It may be worth investigating a number of features, 
notably the treatment of anisotropy and 
multi-layered soils, within the model that have not 
yet been fully tested. , 

An algorithm which varies the time step size should 
be incorporated into SUBFEM. It would reduce the 
time increment during intervals when convergence to 
a solution is expected to be a problem, such as the 
period immediately following changes in boundary 
conditions. The time step size would be increased 
when the numerical system is stable. 

The computing efficiency of SUBFEM could be improved 
by incorporating the substructure method of analysis 
(Elwi and Murray, 1977) into the model. In this 
approach, a hillslope is partitioned into several 


Subunits, each of which is considered to be a large 


finite element. Each subunit in turn is divided into 
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still smaller elements. The larger elements are 
assembled along inter-boundary nodes. Solutions are 
obtained first for the inter-boundary nodes, and 
then for nodes contained within the internal 
Structure of the large elements. The procedure is 
particularly efficient if several large elements are 
identical. 

In many analyses of subsurface flow problems, the 
investigator is interested not in the internal 
Structure of the system, but in conditions at its 
boundaries i.e. at the seepage face or at the flux 
boundaries. For these problems, the Boundary Element 
Method (Brebbia, 1978) may be a more efficient 
technique. It has the advantage that much smaller 
systems of equations are generated, and much less 
data are required to obtain solutions, when compared 
to finite element and finite difference methods. The 
potential application of this method to solving 
subsurface flow problems should be investigated. 
Over the long term, as computer capabilities 
increase and computing costs diminish, it may be 
possible to expand the model to the 
three-dimensional form and, at low cost, use it to 


simulate problems involving entire watersheds. 
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VIII. APPENDICES 


A. Appendix A. Development Of Coordinate Functions Using 
Area Coordinates For Triangular Elements And Integration 
Of Resulting Expressions 
Figure Al is a triangular finite element of area A, 


Meth nodes at i,-j, and k. 


»§ 


Figure Al. Linear triangular finite element and development 
of area coordinates. 


Point P(x,y), which can be any point in the triangle, marks 
the apices of three subtriangles of areas A;, Aj, and A, 
Such that A,+A,+A,=A. The triangular area coordinates of 
P(x,y) are defined as L,=A,/A, Lj=Aj/A, L,=A,/A, and 
consequently L,;+L,+L,=1. Furthermore at node i, L,;=1 and 
L,=L,=0. Similar properties hold at nodes j and k also and 
meet the requirements for defining coordinate functions. 
The area of a triangle can be defined as half the 


Ccross-product of two adjacent sides. Consequently, it can be 


expressed in terms of the triangular coordinates as a 
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determinant of a 3x3 matrix 


1 xX j ¥ i 
2A= {1 xj y; =X j¥e—Xe Vj )— (yey Kit (xe-k yi 
1 X Yr 
or 
2A=X Vu Xe VY) tX1 (Vj -¥u) tyi (x,-x;,) (Ai) 


Similarly, the areas of the subtriangles may be written 


(proceeding in a counterclockwise direction) as 


Lope CY beg Karey eee os: 
2A Sel x AY je) 2A52)|lyxn iyi 3s ) 2Aca |t ey ya 
Pe) y@ eX De Mi Lee 
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2A =(x jp ye Xe Vj )tX (YVR) TY (KE WX; ) (A2.1) 
2A, =( xe yi -Xi Vu) tx Cye yi ty (xix) (A2.2) 
2REE( RY POX iY iP LEAVY OY yA eg Te CA2" 3) 


Since L,;=A,/A, L,=A,;/A, and L,=A,/A, the triangular 


coordinates become 
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These expressions may also be obtained by generating 
coordinate functions using global coordinates. One equation 


which describes a plane over a linear triangular element is 
@,=a,t+b,xtc;\y (A4) 
where a;, b,;, and c; are constants identified with the ith 
coordinate function. Given this definition of a coordinate 
function, together with the prescribed constraints, the 


following set of equations for the ith coordinate function 


is obtained 
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The determinant of the 3x3 matrix is 
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which, as stated earlier, is equal to 2A. The solutions aj, 


b;, and c, are derived from the expressions 
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1 x52 41 
1 X j 0 
1 Xk 0 


ZAC ; = =(x,-x,) 


Thus equation A4 becomes 
@i=CCxXpyeXayj )+(y j-ye) x+(x,-x | yJ/2A 


which is identical to equation A3.1. Similar identities hold 
between ¢; and L;, and between ¢, and L,. 
Having established the relationship between ¢ and L, we 


can proceed to the integration problem 


f,£(¢1,¢),¢,)GA=S, £(L,,L;,L,)dA (A5) 


The integration on the right hand side is easily performed 
if dA is expressed in local coordinates. This can be 


accomplished by using the relation 
aLyaL ,=/J|da 
where |J| is the determinant of the Jacobian matrix 
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and it is a relatively simple matter to compute the six 


types of integrals that appear in the simulation model 


J, Ly da O21, 3,8) (A6.1) 
J, Ly LmdA Cle ik ae, Kee eae (A6.2) 
f, 0b; , dL ;dA (A6.3) 
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B. Appendix B. Program SUBFEM - User's Manual 

The model developed to meet the Ob er tester this 
Study is called the Subsurface Flow Finite Element Model 
(SUBFEM). The computer program, also called SUBFEM, is 
written in FORTRAN, and is run on the Amdahl 470V/8 computer 
at the University of Alberta. Metric units are used 
throughout. The computer code is listed in Appendix C.' 

The first step toward solving the subsurface flow 
problem is to define it. Since the model is limited to 
two-dimensional systems, we are restricted to profiles such 
as those illustrated in figure 10, and must consider the 
factors outlined in Chapter III of this treatise. 

A triangular, finite element mesh is constructed to fit 
the profile. It may be regular (with elements of the same 
size) except near some boundaries, or it may be designed so 
that smaller elements are located in regions where marked 
changes in water content and pressure potential are likely 
to occur, i.e. near the upper boundary and the seepage face. 
Conversely, large elements can be used where such changes 
are expected to be minimal, e.g. in the vicinity of the 
basal, anne boundary. Once the mesh has been constructed, 
Phe nodes are automatically defined since they are located 
at the vertices of the triangles. 

For computational efficiency, the nodes should be 


numbered in the direction of shortest side. In most 


hillslope simulations, this will be in the vertical 


direction, since the profile being Simulated is usually much 
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longer than it is deep. Several examples of finite element 
grids are given in Chapter v. 

Initial and boundary conditions must also be Specified. 
Initially, some assumed or given pressure potential values 
are assigned to each node, thereby establishing initial 
conditions. Three boundary conditions are considered in 
SUBFEM: no-flow,flux, and seepage face. The no-flow boundary 
conditions are applicable to the impermeable boundaries, 
while the flux boundary conditions are usually restricted to 
the upper boundary. Positive fluxes such as rainfall and 
Snowmelt, and negative fluxes such as evaporation, are 
treated in the same manner. Flux values are calculated, and 
assigned to individual flux boundary nodes (Fig. A2). 

It is also necessary to specify the maximum extent of 
the seepage face, even though an unsaturated flow problem is 
being simulated, or though the extent of the seepage face at 
the beginning of the simulation is less than maximum. 

Water withdrawal through transpiration is computed as a 
rate, using the same method employed to compute the fluxes. 
Here, however, the values are assigned to subsurface sink 
nodes, rather than to boundary nodes. 

Information is required on each soil or other material 
that appears in the simulation. First, the Saturated 
hydraulic conductivity as obtained in the field or in the 
laboratory; and secondly, the relation between volumetric 


Soil water content and pressure potential. In the latter 


case, a single-valued, rather than a hysteretic, function is 
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assumed. 

Input data are read in by five subroutines: INPUT1, 
HYCOND, INPUT3, INPUT4, and BOUND. The input data required 
for each of these subroutines, together with the 
identification of variable and parameter names, and the 
specified FORMATS, are described below. 

Input data for cards A to L are stored in a file 
assigned to logical unit #5; they are read in subroutines 


INPUT1, HYCOND, INPUT3, and INPUT4. 


Subroutine INPUT1 reads control data (cards A and B) for the 


system. 


A. Heading Card(1 card) 


HED 
FORMAT (20A4) 
HED sTitle of the problem. 


B. System Card (1 card) 
NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,NS,NC,DELT,NTRANS 


FORMAT(814,F10.5,14) 


NSUB Total number of subsystems. (Insert 1) 
MXND * *Maximum number of inter-boundary nodes 
in any subsystem. (Insert 0) 
NUMNST sTotal number of nodes in the master system. 
MASTR ‘* :Flag for master matrix. (Insert 1) 
MBEL ° Total number of external boundary elements 
attached to master system. (Insert 0) 
IDRY * :If 1, this is a dry run.(Insert 1) 
sIf 0, this is a production run. 
NS *Number of materials (soil types, geological 
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formations, etc.). 

NC :Number of incremented pore classes used to 
compute hydraulic conductivity and specific 
moisture capacity. 

DELT :Time step size. 

NTRANS :If 1, this is a transient problem. 

If 0, this is a steady state problem. 


Subroutine HYCOND reads in cards C to G. 


C. Material Description Card (1 card) 
Sie ol ste, ST3,514,ST5,ST6,ST/ ,STS,STS 
FORMAT (10A8) 


Bit,;ol2,...S19 :Name of material (soil, rock, or till type). 


D. Material Properties Card (1 card) 

INP, TMAX, SCON, SCONR, RESWAT 
FORMAT(15,F10.4,2E15.5,F5.3) 

INP sNumber of input data points for the water 


content/pressure (characteristic) curve - 
limited to 20 points with present format. 


TMAX sMaximum water content (dimensionless) 
SCON :Experimentally obtained saturated 

| hydraulic conductivity (cm/unit time) 
SCONR sRatio of saturated horizontal hydraulic 


conductivity to saturated vertical hydraulic 


conductivity (dimensionless). 
RESWAT *Estimate of residual (immobile) water 
(dimensionless). 


E. Water Properties Card (1 card) 
SURTEN , DENWAT, VISWAT, TEMP, GRAVTY 


FORMAT(F10.6,F10.3,F12.8,2F5. 1) 


SURTEN >Surface tension of water (dynes/cm). 
DENWAT sDensity of water (g/cm*). 
VISWAT sViscosity of water (dyne sec/cm’). 


TEMP sWater temperature (C). 
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GRAVTY Gravitational constant (cm/sec?). 


F, Characteristic Curve (Water Content) Card (1 card) 

PaetAC)),THETA(2),.1. ..+.,THETA(INP) 

FORMAT (20F4.3) 

THETA (1) sData point with lowest water content 
(dimensionless). 


THETA(INP) :Data point with highest water content 
(dimensionless). 


G. Characteristic Curve (Pressure) Card(s) (1 to 3 cards) 

Ber) oP (2); "2 ...,DP(INP) 

FORMAT (8F 10.2) 

DP(1) :Pressure potential that corresponds to THETA(1), 
and which has the highest absolute value (cm). 


DP(INP) :Pressure potential that corresponds to 
THETA(INP), and which has the lowest absolute 


value (cm). 


A set of cards C to G must be included for each of NS 


materials (soils or geologic strata). 


Subroutine INPUT3 reads additional control data for the 


system. 


H. Control Data Card (1 card) 


IS,NNODE,NEL,NSTEP, IMAX,E, SLOPE, NFIBN,NIBNS ,NSEEP 


FORMAT (514,2E12.5,314) 


IS :Subsystem number. 

NNODE sNumber of nodes 1n the system. 
NEL sNumber of elements in the system. 
NSTEP sNumber of time steps. 


IMAX sMaximum number of iterations per time step. 
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E :Specified difference criterion for 
iterations (cm). 

SLOPE sHillside slope (degrees). 

NF IBN * Local number of first inter-boundary node. 

NIBNS * :Number of inter-boundary nodes. 

NSEEP sNumber of seepage face nodes. 


Subroutine INPUT4 reads and writes the nodal and element 


data, together with initial pressure potential. 


I. Nodal Geometry Cards (i card per node) 
N,X(N) ,Z(N) 


FORMAT (I16,2F12.2) 


N s:Node number. 

X(N) :X coordinate of node. 
Z(N) :Z coordinate of node. 
Notes: 


1)Nodal data cards must be entered in numerical order. If 
nodes are omitted then the program generates coordinates for 
intermediate nodal points. 

2)For constant pressure potential nodes, pressure 
potential(w) values are known, and the matrices are 
partitioned so that equations for these nodes are 
eliminated. Boundary condition information related to these 


nodes is transferred to the known matrix. 


J. Initial Pressure Potential Card (1 or more cards) 
H(1),H(2),... «..,H(NNODE) 


FORMAT (10F8.2) 


H(1) sInitial pressure potential at node I 


K. Material Identity Card (1 or more cards) 
MATL(1),MATL(2),... .. «MATL (NNODE ) 


FORMAT (4012) 


on” 


to deena 3  sapaannae be 
iy ,. : a 


\eion yssbaved= geet GR to secre 


*~ 


- or a 

rasnels bas (shor ade Wetiaw Sos ehkey STOeK eam 
Ww ~ ‘“s “F . ee 

‘stsnetog S4¥URBese Laisinl atiw yeddepe 


_ 


(ston) taq Br62 |) ahasd 3180080 | 
4 
¥ 


. one : y 
(S.SteG ah 


.tsdayn ebous. = 
;a8on io stshihioes Be ie 
~f6on 30 2tsnibugoas | 
ey | 

: | 7 ol 


ii .%9090 \teosasmrs aE bBseistns sd Jigga atis>: sat fal 
103 zasanib 60d sedbtonap mexporg edd neds bedding: 
|. .etniog Labon stsibem 
| Saidnesoq-e1vess aaaser 09.5 
S ,QWOGSs 938 aay sit 
“BOR eaoiteups t9d2 o@ 
Ke Mattihsop ya 
| a aris ea beste. 


pxvReSIo za 
216 gepisiem 
, ste 8abon At 
s2ed+ o3 bstslet fod 


TX 


on 


(ebigo 08 ahiys cis Asitzas309 otiaaees ts int 
ire mi ee 
or Be aba (BQOUW) Hi. ws areal 
my Piet 5 (senor 


tate $8 ‘staan aayeesag ratdiand ‘aie . 


Sieh 
yam fai 


CS, 


MATL (I) :Code number for material (soil, till, sandstone, 


etc.) at node I. 


L. Element Data Cards (1 card per element) 


M,NP(1,M),NP(2,M) ,NP(3,M) 


FORMAT (416) 
M :Element number. 
NP(I,M) sNodal point numbers at the three corners of 


the triangular element in counterclockwise order. 
Notes: 
Element cards must be in order. If element numbers are 


omitted, the program generates the intermediate element 
data. 


Data for cards M to V are stored in a file assigned to 


logical unit #7; they are read in subroutine BOUND. 


M. New Boundary Conditions Card (1 card) 

NEWBC 

FORMAT (13) 

NEWBC eNumber of times boundary conditions change. 


N. Boundary Conditions, Time-of-change Card (1 card) 
KSTNO(1),KSTNO(2),... ««+,KSTNO(NEWBC ) 


FORMAT (2513) 


hichetinse 
KSTNO (1 sNumber of time step during w 
Bee change in boundary conditions occurs. 


h last 
sNumber of time step during whic 
ui change in boundary conditions occurs. 


O. Seepage Face Nodes Card (1 card) 
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NSPF(1),NSPF(2),... ...,NSPF(NSEEP) 


FORMAT (2513) 


NSPF (1) :Nodal number of first seepage face node. 
NSPF(NSEEP) :Nodal number of last seepage face node. 


Note: 


If there is no outflow at the beginning of a simulation, 
this card is omitted. 


P. Specialized Nodes Identification Card (1 card) 


NEUM,NFLUX,NSINK,NDIR 

FORMAT (415) 

NEUM :Total number of no-flow boundary nodes. 

NF LUX :Total number of flux boundary nodes. 

NSINK :Total number of sinks. 

NDIR :Total number of Dirichlet (constant pressure 


potential) boundary nodes. 


Q. Boundary Conditions and Sinks Card#1 (1 to 3 cards) 
Benin) BYCD( 2) ,.\9. .« > ,BYCDV20) 


FORMAT (8E10.3) 


BYCD(1),... ..-,BYCD(5):Flux values to be assigned to flux 


boundary nodes’ 
BYCD(6),... ...,BYCD(10):Withdrawal rates to be assigned to 


Sink -nodes’. 
BYCD(11),... ...,BYCD(20):Pressure potential values to be 
assigned to constant pressure potential (Dirichlet) nodes’. 


Note: 


One BYCD value may be assigned to several boundary (or sink) 
nodes. 


R. Boundary Conditions and Sinks Card#2 (1 card) 


NECCNA NEG (2), JPL PES, NBC (20° 


7Identified on cards T through V. 
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FORMAT (2014) 


NBC(KZ):Cumulative number of nodes to which values BYCD(Kz) 
must be assigned. 

NBC(1),... ...,NBC(5):Flux nodes. 

NBC(6),... ...,NBC(10):Sink nodes. 

NBC(11),... ...,NBC(20):Constant pressure potential 
(Dirichlet) nodes. 


Example; If BYCD(1)=20.0, BYCD(2)=40.0, BYCD(3)=50.0, 

and NBC(1)=5, NBC(2)=8, and NBC(3)=18, 

then the value 20.0 is assigned to the first 5 

flux nodes, the value 40.0 to the next 3 (8 minus 5) . 
flux nodes, and the value 50.0 to the next 10 (18 minus 8) 
flux nodes. 


S. No-flow Boundary Nodes - Identification Cards (1 or more 
cards) 

NBOUND(1),NBOUND(2),... ...,NBOUND(NEUM) 

FORMAT (2014) 


NBOUND( 1) :Number of first no-flow boundary node. 
NBOUND(NEUM) :Number of last no-flow boundary node. 


T. Flux Nodes - Identification Cards (1 or more cards) 
NBOUND(1),NBOUND(2),... ..+,NBOUND(NFLUX) 


FORMAT (2014) 


NBOUND ( 1 ) sNumber of first flux boundary node. 
NBOUND(NFLUX) :Number of last flux boundary node. 


U. Sink nodes - Identification Cards (1 or more cards) 
NBOUND(1),NBOUND(2),... . » ,NBOUND (NSINK) 


FORMAT (2014) 


NBOUND ( 1 ) sNumber of first sink node. 
NBOUND(NSINK) :Number of last sink node. 
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V. Dirichlet (Constant Pressure Potential) Nodes - 
Identification Cards (1 or more cards) 
NBOUND(1),NBOUND(2),... ...,NBOUND(NDIR) 

FORMAT (2014) 


NBOUND ( 1) :Number of first Dirichlet boundary node. 
NBOUND(NDIR) :Number of last Dirichlet boundary node. 


An example of a complete set of input data for SUBFEM is 


listed in Appendix D. 


Subroutines HYCOND and TABLOK are used to compute 
hydraulic conductivity and specific moisture capacity, as a 
function of pressure potential, for each material in the 
system. The hydraulic conductivity is calculated using a 
program developed by Green and Corey (1971). This 
computational procedure is based on knowledge of saturated 
hydraulic conductivity, and on the relation between 
unsaturated hydraulic conductivity and pore-size 
distribution. TABLOK is adapted from the subroutine of the 
Same name used in the TEHM model (Huff et al, 1977). 

Subroutine GMATX solves the integral equations and 
forms the element matrices. It subsequently assembles 
contributions from each element matrix into the global 
matrices [A] and [B]/At. Finally it produces the terms 
[A]+[B]/At and [B]{P}/At. 

The boundary conditions are handled by subroutine 
BOUND. It assigns values to the flux and sink nodes, and 


partitions from the global matrix the equations for the 
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constant pressure potential nodes. 

Subroutine SKYLIN locates the first non-zero element in 
each column of the global matrix,and eliminates all the zero 
elements above that element. This procedure improves storage 
efficiency. SKYLIN also keeps account of column lengths and 
the locations of the diagonal elements in the modified 
matrix array. 

Decomposition of the global matrix using Cholesky's 
method is instituted in subroutine DECOMP. Backsubstitution 
to obtain the solutions to the specified problem is carried 
out in subroutine BKSB1. 

Fluxes are computed in subroutine FLUX. This subroutine 
also determines the extent of the seepage face and the 
ouct low from it. 

Subroutine ITERAT controls the iteration process and 
prints out the results at the end of each time step. Results 
include the fluxes, the pressure potential field, the number 
of iterations completed during the time step, and a measure 
of convergence. 

Subroutine XPLOT stores into a file (assigned to 
logical unit #8) for plotting purposes, information 
pertaining to the final time step. Data stored for each node 
include: the coordinates, pressure potential, total 
potential, volumetric water content, and the hydraulic 
conductivity. 

Plotting is done on the AED 512 colour graphics 


terminal, using the IGPLSD integrated graphics subroutine. 
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This arrangement permits three-dimensional display of data, 
and facilitates data manipulation for plotting purposes. 
Hardcopy is obtained from the University of Alberta's 


CalComp 925/1036 plotter. 
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C. Appendix C. Program SUBFEM - Computer Code 


SUBFEM 
1 C SUBSYSTEM ANALYSIS OF FLOW THROUGH POROUS MEDIA, 
BU gps tds don kis lovaehans cance eye ee Nees 
Linen Ree eRe en Dees Pekakeeeee 
4 IMPLICIT REAL*8(A-H,0-Z) 
5 REAL*8 NAMES 
€ Cc 
7 COMMON /MASTIA/ NNN( 1000) 
8 COMMON /MASTRA/ BBB( 1000) 
g COMMON /SUBIA/  MMM(3500) 
10 COMMON /SUBRA/ AAA(85000) 
11 COMMON /SUBSR/ CCC(4000) 
12 COMMON /DIMCOM/ L1,L2,L3,L4,L5,MXDM,NAMES(5,20),IPT(5,21), 
13 *ICOM(5) 
14 Cc 
15 ; ICOM(1)=1000 
16 ICOM(2)=1000 
17 ICOM(3) =85000 
18 ICOM( 4) =3500 
19 ICOM(5)=4000 
20 c 
21 CALL MAINMG 
22 c 
23 . END 
24 C 
25 Cc 
26 Cc 
OH, Cc 
28 c 
29 SUBROUTINE MAINMG 
30 € THIS SUBROUTINE MANAGES THE READING, THE FORMULATION, 
31 Cc AND THE SOLUTION OF THE PROBLEM. 
a RS PELOTON SE LAE OPE E EEE BE REESE ye Gr 
34 IMPLICIT REAL*8(A-H,0-Z) 
35 REAL*8 NAMES,NAME,KZZ,KRATIO 
36 C 
37 COMMON /BORDER/ NEUM,NDIR,KFLAG, NBOUND(60) 
38 COMMON /HILL/ SLOPE,COS1,COS2,SIN1,SIN2 
39 COMMON /MASTIV/ NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,IN,10,DELT,NTRANS 
40 COMMON /MASTIA/ NNN( 1000) 
41 COMMON /MASTRA/ BBB( 1000) 
42 COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 
43 COMMON /SUBIA/ MMM(3500) 
44 COMMON /SUBRA/ AAA(85000) 
45 COMMON /SUBSR/ CCC(4000) 
46 COMMON/HOOKUP/NS ,NC, RESWAT 
47 COMMON /DIMCOM/ LA1,LA2,LA3,LA4,LA5,MXDIM,NAMES(S, 20), 
48 *IPT(5,21),1ICOM(5) 
49 C 
50 CH AH RANT OR IR RRR He PRC RAR AOE AE AE OE NE eT Sets 
51 IN=5 
52 10=6 
53 REWIND 1 
54 ND=1 : 
55 NST=1 
56 ITTER=O 
57 KFLAG=O 


58 NEW=O 
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READ SYSTEM CONTROL VARIABLES 


CALL INPUT 14 


CALCULATE HYDRAULIC CONDUCTIVITY AND SPECIFIC MOISTURE CAPACITY 
AS A FUNCTION OF PRESSURE HEAD 


I 1=ISPAC(6HWATCON,NS*NC, 1) 
I2=ISPAC(SHSLPOT,NS*NC, 1) 

I3=ISPAC(SHSLCON,NS*NC, 1) 

14=ISPAC(5HSLCAP,NS*NC, 1) 

15=ISPAC(G6HKRATIO,NS, 1) 

CALL HYCOND(888(11),88B(12),888(13),BBB(14),B8BB(I5)) 


READ CONTROL DATA 
30 CALL INPUT3(NSTEP, IMAX,E) 


READ SYSTEM GEOMETRY, ELEMENT CONNECTIVITY AND PROPERTIES 


K1=ISPAC(1HX,NNODE, 3) 
K2=ISPAC(1HZ,NNODE, 3) 
.K3=ISPAC(1HH,NNODE,3) 

K4=ISPAC(3HPSI ,NNODE,3) 
KS5=ISPAC(4HPSIP,NNODE, 3) 
K6=ISPAC(SHTHETA,NNODE, 3) 
K7=ISPAC(3HKZZ,NNODE, 3) 
K8=ISPAC(2HCA,NNODE, 3) 
K9=ISPAC(4HGMAT , NNODE*(NNODE+1)/2,3) 
K 10=ISPAC(4HBMAT , NNODE*(NNODE+1)/2,3) 
K11=ISPAC(3HRHS,NNODE, 3) 
K12=ISPAC(4HOQRHS ,.NNODE , 3) 
K13=ISPAC(4HSINT,NNODE, 3) 
K14=ISPAC(SHSINT2,NNODE,3) 
K15=ISPAC(4HDIFF,NNODE,3) 
K16=ISPAC(1HQ,NNODE, 3) 
K17=ISPAC(4HTPOT,NNODE,3) 
L41=ISPAC(2HNP,NEL*3, 4) 
L2=ISPAC(4HMATL ,NNODE, 4) 
L3=ISPAC(2HLD,NNODE, 4) 
L4=ISPAC(3HLDQ,NNODE ,4) 
L5=ISPAC(4HKODE ,NNODE, 4) 
LG=ISPAC(4HNSPF ,NNODE,4) 
L7=ISPAC(SHKOLHT ,NNODE, 4) 

CALL INPUT4(AAA(K1),.AAA(K2) ,AAA(K3) , AAA(K4), 


*MMM(L1),MMM(L2)) 


READ PRESSURE POTENTIALS, WATER CONTENTS, CONDUCTIVITIES, 


AND CAPACITIES FROM TABLES 
SoeCALE TABLOK(AAA(K3) , AAA(K4) , AAA(KG&) , AAA(K7) , AAA(K8), 


*BBB(1I1),8BB(I2),88B8(13),8BB(14),MMM(L2)) 


IF(NST.LE.NSTEP) GO TO 150 
CALL XPLOT(AAA(K1) , AAA(K2) , AAA(K3) , AAA(K6) , AAA(K7), 


*AAA(K1i7)) 


FORM SYSTEM MATRIX FROM ELEMENT MATRICES 
150 CALL GMATX(AAA(K1),AAA(K2), AAA(K3) , AAA(K7) , AAA(K8) , AAA(K9), 


*AAA(K10) ,AAA(K11),AAA(K12) ,MMM(L1),MMM(L2),MMM(L3),MMM(L4), 
*BBB(I5)) 
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DEFINE BOUNDARY CONDITIONS 
CALL BOUND(AAA(K9) ,AAA(K11),AAA(K12),AAA(K13), 
*AAA(K14),MMM(L3),MMM(LS),MMM(L6),NSTEP, NST ,NEW) 


FORM DIAGONAL ADDRESSING ARRAY AND 
COMPUTE COLUMN HEIGHTS 
CALL SKYLIN(AAA(KS) ,MMM(L3),MMM(L7)) 


CHOLESKY DECOMPOSITION 
CALL DECOMP(AAA(KQ) ,AAA(K13),MMM(L7), MMM(L3) ) 
CALL BKSB1(AAA(K9),AAA(K13), MMM(L3))_ 


CARRY OUT ITERATIONS 
CALL ITERAT(AAA(K1) ,AAA(K2),AAA(K3) , AAA(K4) , AAA(K5) , AAA(K10), 
*AAA(K11),AAA(K12), AAA(K13), AAA(K15)., AAA(K16), AAA(KI7), MMM(L4), 
*MMM(L5),MMM(L6),IMAX,E,NST.NSTEP) 
GO TO 85 
100 RETURN 
END 


SUBROUTINE BKSB1(GMAT,SINT,LD) 


BACKSUBSTITUTION ALGORITHM: SOLVES THE PROBLEM 
D*L(T)*U=R FOR MASTER SYSTEM (EITHER BOUNDARY 
NODES MATRIX, OR UNSUBSTRUCTURED SYSTEM). 


Se KK oe eke ok ok ok ok ak ok ok ok ot ok ok ok ke ke ake ak ake oie ok ake ke aie aK ok ak ake ac ok oe ke ake ke ake ke ke ke ok ke ke oi ake fe i oa a ok fe ok ae 


IMPLICIT REAL*8(A-H,0-Z) 
COMMON /MASTIV/ NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,IN,IO,OELT,NTRANS 


COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 


DIMENSION GMAT(1),SINT(1),LD(1) 


COMPUTE D(-1)*R 


IF(MASTR) 250,250, 100 
100 INTNOD=NTEMP 


WRITE(6,500) 
500 FORMAT(//,10X, “POTENTIAL FIELD AND SURFACE INTEGRALS FOR 


*MASTER SYSTEM’ ,//) 
DO 200 JG=1, INTNOD 


L=LD(uG) 
SINT (JG) = =SINT(JUG)/GMAT(L) 


200 CONTINUE 
GO TO 50 


BACKSUBSTITUTION U=eGh) (1) *R 


250 INTNOD=NFIBN- 1 


WRITE(6,600) 
GOO FORMAT(//,10X,’GMAT AND PRESSURE HEAD FIELD’ ,//) 


50 DO 300 JG=2,INTNOD 
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THIS SUBROUTINE DEFINES THE BOUNDARY CONDITIONS FOR THE 
SYSTEM , ASSIGNS VALUES TO THE SURFACE INTEGRALS, AND 
DELETES EQUATIONS FOR THE CONSTANT HEAD NODES 
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IF BOUNDARY CONDITIONS REMAIN UNCHANGED,REASSIGN ORIGINAL VALUES 


N=INTNOD-uJG+2 
KL=LD(N-1)+14 

KU=LD(N) -1 

IF(KL.GT.KU) GO TO 300 
K=KL-LD(N)+N 

DO 350 KK=KL,KU 
SINT(K)=SINT(K)-GMAT (KK) *SINT(N) 
K=K+1 

CONTINUE 

IALL=LD(INTNOGO) 

WRITE(6,700) (GMAT(I),1I=1,IALL) 
FORMAT(10E 12.4) 

WRITE(6,700) (SINT(I),1I=1, INTNOD) 
RETURN 

END 


BLOCK DATA 


REAL*8 NAME ,NAMES 


COMMON /DIMCOM/ LAST1,LAST2,LAST3,LAST4,LASTS,MAXDIM, 
*NAMES(5,20),IPT(5,21),ICOM(5) 


DATA LAST1,LAST2,LAST3,LAST4,LAST5,MAXDIM, IPT(1,1),IPT(2,1), 
zTPT( 3, 1), 1PT(4, 1), 1PTCS, 1)/5*0, 20,5717 


END 


SUBROUTINE BOUND(GMAT,RHS,QRHS,SINT,SINT2,LD,KODE, 
*NSPF,NSTEP,NST,NEW) 


IMPLICIT REAL*8(A-H,0-Z) 


COMMON /BORDER/ NEUM,NDIR,KFLAG, NBOUND(60) 
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COMMON /MASTIV/ NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,IN,IO,DELT,NTRANS 


COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 


DIMENSION GMAT(1),SINT(1),KODE(1),L0(1),NSPF(1),RHS(1), 
*SINT2(1),BYCD(20) ,NBC(20) , QRHS(1),KSTNO(5O) 


IALL=NNODE*(NNODE+1)/2 
NTEMP=NNODE 

WRITE(10, 1300) 

WRITE(IO, 200) (GMAT(I),1=1,1ALL) 


FORMAT(//,5X, ‘THE GLOBAL (G) MATRIX IS:’//) 
IF(NST.EQ.1.AND.ITTER.EQ.O) GO TO 65 


TO SURFACE INTEGRAL ARRAY, SINT(I) 
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239 Cc 
240 DO 150 KL=1,NNODE 

244 150 SINT(KL)=SINT2(KL) 

242 IF(ITTER.EQ.0O) GO To 75 

243 IF(NDIR.EQ.0) GO TO 700 

244 K=4 

245 NEND=NDIR 

246 GO TO 95 

247 Cc 

248 Cc READ AND WRITE INITIAL CONDITIONS (NUMBER OF NO-FLOW,FLUX, 
249 Cc SINK, CONSTANT-HEAD,SEEPAGE NODES,ETC) 
250 C 

251 65 NPLUS=1 

252 READ(7,750) NEWBC 

253 READ(7,22) (KSTNO(KBC),KBC=1,NEWBC) 
254 WRITE(IO,22) (KSTNO(KBC),KBC=1,NEWBC) 
255 IF(NSEEP.EQ.0) GO TO 4000 

256 READ(7,22) (NSPF(L),L=1,NSEEP) 

257 WRITE(1I0,22) (NSPF(L),L=1,NSEEP) 

258 - 4000 READ(7,50) NEUM,NFLUX,NSINK,NDIR 
259 WRITE(I0O,50) NEUM,NFLUX,NSINK,NDIR 
260 DO 90 I=1,NNODE 

261 80 KODE(I1)=0 

262 K=1 

263 NEND=NEUM 

264 KA=1 

265 KB=6 

266 GO TO 85 

267 Cc 

268 C DETERMINE IF NEW BOUNDARY VALUES ARE TO BE READ IN 
269 G 

270 75 NEW=NST-KSTNO(NPLUS) 

271 IF(NEW) 220,221,221 

Mo 220 IF(NDIR.EQ.0O) GO TO 700 

273 K=4 

274 NEND=NDIR 

275 GO TO 95 

276 221 NPLUS=NPLUS+1 

2a Cc 

278 Cc READ AND WRITE INITIAL OR NEW BOUNDARY VALUES 
279 Cc 

280 85 DO 905 I=1,NNODE 

281 805 SINT(I)=0.0 

282 READ(7,33) (BYCD(KY),KY=1, 20) 

283 WRITE(1I0,33) (BYCD(KY),KY=1, 20) 

284 READ(7, 100) (NBC(KZ),KZ=1,20) 

285 WRITE(I0, 100) (NBC(KZ),KZ=1,20) 

286 IF(NST.EQ.1) GO TO 400 

287 GO TO 810 

288 400 READ(7,100) (NBOUND(I),1=1,NEND) 

289 WRITE(IO, 100) (NBOUND(I),1=1,NEND) 
290 IF (NST.EQ.1) GO TO 95 

291 G 

292 CG ASSIGN NEW BOUNDARY VALUES TO BOUNDARY NODES 
293 (3 

294 810 IA=O 

295 IB=5 

296 KA=1 

297 KB=6 


298 DO 3500 J=1,NNODE 
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IF(KODE(JU)-11) 3500,3600,3700 
3600 IA=IA+14 
IF(IA-NBC(KA)) 1,1,2 
2 KA=KA+1 
1 SINT(J)=BYCD(KA) 
GO TO 3500 
3700 IB=IB+1 
IF(IB-NBC(KB)) 5.5,6 
6 KB=KB+1 
5 SINT(J)=BYCD(KB) 
3500 CONTINUE 
K=4 
NEND=NDIR 


ASSIGN INITIAL VALUES TO BOUNDARY NODES 


95 KC=11 
DO 600 I=1,NEND 
J=NBOUND (I) 
GO TO (10,20,30,40),K 
10 KODE(J)=1 
IF(I.NE.NEND) GO TO 600 
NEND=NFLUX 
K=2 
IF(NEND) 15,15,650 
20 IF(I-NBC(KA)) 1900, 1900, 2000 
2000 KA=KA+14 
1900 SINT(J)=BYCD(KA) 
KODE(J)=11 
IF(I.NE.NEND) GO TO 600 
6G0_TO 15 
30 IF(I-NBC(KB)) 2100,2100,2200 
2200 KB=KB+1 
2100 SINT(J)=BYCD(KB) 
KODE (J)=12 
IF(I.NE.NEND) GO TO 600 
GO TO 80 
415 NEND=NSINK 
K=3 
IF(NEND) 80,80,650 
80 NEND=NDIR 
K=4 
IF(NEND) 140, 140,650 
40 IF(NEW.LT.O) GO TO 2400 
IF(ITTER.GT.O.OR.I.GT.1) GO TO 2400 
140 DO 1000 MN=1,NTEMP 
1000 SINT2(MN)=SINT(MN) 
IF(NDIR.EQ.O) GO TO 700 
2400 IF(I-NBC(KC)) 1400, 1400, 1800 
1800 KC=KC+1 
1400 KODE(J)=10 


250 


CONDENSE MATRIX BY ELIMINATING EQUATIONS FOR CONSTANT HEAD NODES 


JCUR=J-I+1 
JJ=JCUR- 1 
IF(JUJ) 500,500,800 


TRANSFER CONSTANT HEAD TERMS TO RIGHT HAND SIDE OF EQUATIONS 
AND DELETE EQUATIONS FOR CONSTANT HEAD NODES 
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359 Cc 
360 800 DO 60 KL=1,Ju 

361 KK=LD(JUCUR) -JCUR+KL 

362 60 RHS(KL)=RHS(KL)-BYCD(KC) *GMAT (KK) 

363 500 NN=JCUR+14 

364 DO 300 LL=NN,NTEMP 

365 MM=LD(LL)-LL+UCUR 

366 RHS(LL)=RHS(LL)-BYCD(KC) *GMAT (MM) 

367 RHS(LL-1)=RHS(LL) 

368 300 SINT(LL-1)=SINT(LL) 

369 KIN=LD(JCUR)+1 

370 IALL=NTEMP*(NTEMP+1)/2 

371 DO 70 II=KIN,IALL 

372 70 GMAT(II-JUCUR)=GMAT(II) 

373 NCOLS=NTEMP-UCUR 

374 DO 110 KI=1,NCOLS 

375 LIN=LD(JCUR+KI-1)+1 

376 LA=LD(UCUR+KI)-1 

377 DO 120 II=LIN,LA 

378 120 GMAT(II-KI)=GMAT(II) 

379 110 CONTINUE 

380 NTEMP=NTEMP- 1 

381 IF(I.EQ.NEND) GO TO 700 

382 600 CONTINUE 

383 650 CONTINUE 

384 GO TO 400 

385 700 IALL=NTEMP*(NTEMP+1)/2 

386 C1100 WRITE(1I0,5000) 

387 5000 FORMAT(//,5X,’THE SINT TERMS (BO) ARE:’,//) 
388 Cc WRITE(IO,200) (SINT(IA),IA=1,NTEMP) 

389 (6 WRITE(10,5500) 

390 5500 FORMAT(//,5X,’THE RHS TERMS (BO) ARE:’,//) 
391 C WRITE(I0,200) (RHS(IB),IB=1,NTEMP) 

392 1100 DO 900 MN=1,NTEMP 

393 800 SINT(MN)=SINT(MN)+RHS(MN) 

394 C WRITE(IO, 1500) 

395 14500 FORMAT(//,5X,’THE CONDENSED GLOBAL (G) MATRIX IS:’,//) 
396 (es WRITE(IO,200) (GMAT(I),1=1,IALL) 

397 Cc WRITE(I0O, 1600) 

398 Cc WRITE(10,200) (RHS(I),1=1,NTEMP) 

399 1600 FORMAT(//,5X,’THE CONDENSED RIGHT HAND SIDE IS:’,//) 
400 Cc WRITE(IO,200) (SINT(I),I=1,NTEMP) 

401 Cc WRITE(1I0, 100) (KODE(I),I=1,NTEMP) 

402 Cc WRITE( 10,2500) ; 

403 2500 FORMAT(///,’THE BOUNDARY CONDITIONS FOR THE NEXT ITERATION ARE:’,/ 
404 1//) 

405 (ie WRITE(I0,200) (SINT2(MN) ,MN=1,NNODE ) 

406 RETURN 

407 C 

408 C FORMAT STATEMENTS 

409 Cc 

410 22 FORMAT(2513) 

411 33 FORMAT(8E10.3) 

aa 50 FORMAT(415) 

413 1400 FORMAT(2014) 

414 200 FORMAT(10E12.4) 

415 750 FORMAT(I3) 

416 END 

447 C 


418 Cc 
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SUBROUTINE DECOMP(GMAT,SINT,KOLHT,LD) 


Cc 
Cc 
Cc CHOLESKY DECOMPOSITION-~ALGORITHM, TO COMPUTE G(I,UJ),L(T), 
Cc AND D IN THE L(T) *D DECOMPOSITION,IT ALSO COMPUTES THE 

Cc RI* VECTOR AND PROVIDES FOR SKYLINE CONFIGURATION 

Cc 

Cc 

C 
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IMPLICIT REAL*8(A-H,0-Z) 
@ 

COMMON /MASTIV/ NSUB,MXND,NUMNST,MASTR,MBEL, IDRY,IN,10,DELT,NTRANS 
C 

COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP, NTEMP 
C 

DIMENSION GMAT(1),SINT(1),KOLHT(1),LD(1) 

ITALL=LD(NTEMP) 
Cc WRITE(IO,300) (GMAT(i),I=1,IALL) 
é WRITE(I0,300) (SINT(I),I=1,NTEMP) 

IF (KOLHT(2).£Q.1) GO TO 50 

GMAT (2) =GMAT(2)/GMAT(1) 

SINT(2)=SINT(2)-GMAT(2)*SINT(1) 

GMAT (3) =GMAT(3)-(GMAT(1)*GMAT(2) **2) 


FORM ELEMENTS OF GMATRIX 


aaa 


50 IF(MASTR) 100, 100,200 
200 INTNOD=NTEMP 
Cc NNODE =NUMNST 

GO TO 150 

100 INTNOD=NFIBN-1 

150 DO 400 JG=3, INTNOD 
IF (KOLHT(JG).EQ.1) GO TO 400 
MJ=JG-KOLHT (JUG) +14 
IF(MU+1.GT.JUG-1) GO TO 550 
MJ1=MJ+1 
JG 1=JG- 1 
DO 500 IG=MJi,JG1 
MI=IG-KOLHT(IG)+1 
IF(MJ-MI) 220,220,230 

230 MI=MdJ 

220 SUMPR=0.0O 
IGi=IG-1 
DO 600 IGS=MI,IG1 
M1=LD(IG)-IG+IGS 
M2=LD(UG)-JG+IGS 
SUMPR=SUMPR+GMAT (M1) *GMAT (M2) 

600 CONTINUE 
MM=LD(UG)-JUG+IG 
GMAT (MM) =GMAT (MM) -SUMPR 

500 CONTINUE 


FORM ELEMENTS FOR L(T) MATRIX AND COMPUTE RI* VECTOR 


aaa 


550 SUMPR=0.0 
JG1=JG- 1 
DO 650 IG=MUJ,JG1 
L=LD( 1G) 
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MM=LD(JG)-UG+IG 
GMAT (MM) =GMAT (MM) /GMAT(L) 

650 SUMPR=SUMPR+GMAT (MM) *SINT(IG) 
SINT(JUG)=SINT(JUG)-SUMPR 


Cc 
Cc FORM DIAGONAL ELEMENTS 
Gc 

SUMPR=0.0 

JG1=JG-1 

DO 700 IGS=MJ,UG1 

LL=LDICIGS) 


MD=LD( JUG) -JUG+IGS 
SUMPR=SUMPR+GMAT (LL) *GMAT (MD) **2 
700 CONTINUE 
MD=LD( UG) 
GMAT (MD) =GMAT(MD)-SUMPR 
400 CONTINUE 


Cc WRITE(IO,300) (SINT(I),1I=1,NTEMP) 
IALL=LO(NTEMP) 

c WRITE(IO,300) (GMAT(I),I=1,IALL) 
RETURN 


Cc 
Cc FORMAT STATEMENTS 
300 FORMAT(10E12.4) 
END 


SUBROUTINE FLUX(X,Z,PSI,BMAT,QRHS,SINT,Q,TPOT,LDQ,KODE,NSPF ) 


Cc 

Cc THIS SUBROUTINE COMPUTES FLUXES, DETERMINES THE EXTENT OF THE 
Cc SEEPAGE FACE, AND THE OUTFLOW FROM IT 
Cc 
C 
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IMPLICIT REAL*8(A-H,0O-Z) 

COMMON /BORDER/ NEUM,NDIR,KFLAG,NBOUND(60) 

COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 
DIMENSION PSI(1),BMAT(1),SINT(1),Q(1),L0Q(1),KODE(1), 
*x(1),Z(1),TPOT(1),QRHS(1),NSPF(1) 


Cc 
Cc INITIALIZE @ ARRAY AND COMPUTE HYDRAULIC HEAD 
ce; 
KFLAG=O 
DO 110 I=1,NNODE 
Cc TPOT(1)=Z(1)+PSI(I) 


110 Q(1)=0.0 
WRITE(8, 850) NNODE 
WRITE(6, 300) 


COMPUTE FLUXES 


agaaNaANAaADA 


DO 1600 I=1,NNODE 
KOLT=LOQ(1)-I 

KDIAG=LDQ(1I) 

DO 500 J=1,NNODE 

KOUNT =KOL T+J 
IF(KOUNT.GT.KDIAG) GO TO 650 
Q(1)=Q(1)+BMAT(KOUNT) *PSI(J) 
GO TO 500 
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650 


500 


1600 


2000 


254 


JM=J-1 
KROW=LDQ(UM)+I 
Q(1)=Q(1)+BMAT(KROW) *PSI (J) 

CONT INUE 

Q(1I)=Q(1)-ORHS(1) 

WRITE(8,950) I,X(I),Z(1),PSI(1I),TPOT(I) 
CONT INUE 

WRITE(6,100) (PSI(I),1I=1,NNODE) 

WRITE(6, 200) 

WRITE(6,100) (Q(I),I=1,NNODE) 

WRITE(6, 2000) 

FORMAT(//,5X,’THE QRHS TERMS (FL) ARE:’,//) 
WRITE(6,100) (QRHS(I),I1=1,NNODE) 
WRITE(6,50) 


IF(NSEEP.EQ.0O) GO TO 1000 
DO 600 MI=1,NSEEP 
T=NSPF(MI) 

IF(KFLAG) 900,900,800 


IDENTIFY CONSTANT HEAD BOUNDARY NODES 


9300 


IF(KODE(I).EQ.1) GO TO 750 
IF(Q(1I).LE.0.0000001) GO TO 600 


MAKE NODE A "NO-FLOW" BOUNDARY NODE 


800 


Uso 


IF(KODE(I).EQ@.1) GO TO 750 

Q(1)=0.0 

KODE(I)=1 

NEUM=NEUM+ 14 

NDIR=NDIR- 1 

KFLAG=1 

WRITE(6,150) I,KODE(I),PSI(I),Q(1),NEUM,NDIR 
GO TO 600 

TFCPSI(1).LT.O.Q) GO TO GOO 


MAKE NODE A CONSTANT HEAD BOUNDARY NODE 


700 


600 


PSI(1I)=0.0 

KODE(1I)=10 

NDIR=NDIR+1 

NEUM=NEUM- 4 

WRITE(6, 150) I,KODE(I),PSI(1),Q(1),NEUM,NDIR 


CONTINUE 


REDEFINE SEEPAGE FACE 


250 
1000 


JJ=0 

DO 250 J=1,NNODE 
IF(KODE(JU).NE.10) GO TO 250 
JJ=Ju+1 

NBOUND (JJ) =J 

CONTINUE 

CONTINUE 

RETURN 


FORMAT STATEMENTS 
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299 
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255 


5O FORMAT(//,10X,’NODE’,10X, ‘KODE’, 10X, ‘PRESSURE’, 14X, ‘FLUX’, //) 
100 FORMAT(10E12.4) 

150 FORMAT(10X,14,10X,14,8X,E£12.5,8X,E£12.5,215) 
200 FORMAT(//,5X,’THE COMPUTED FLUXES ARE’ ,//) 
300 FORMAT(//,5X,’THE PRESSURE HEAD FIELD IS:’,//) 
850 FORMAT(14) 
950 FORMAT(14,4F12.1) 

END 


aan 


SUBROUTINE GMATX(X,Z,H,KZZ,CA,GMAT,BMAT,RHS,QRHS,NP,MATL, 
*LD,LDQ,KRATIO) 


THIS ROUTINE FORMS THE ELEMENT MATRICES FOR EACH TRIANGULAR ELEMENT, 
SOLVES THE INTEGRAL EQUATIONS,AND ASSEMBLES THE CONTRIBUTIONS 
FROM EACH ELEMENT COEFFICIENT MATRIX TO FORM THE GLOBAL MATRICES. 
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IMPLICIT REAL*8 (A-H,O-Z) 

REAL*8 KRATIO,KZZ 

COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 
COMMON /MASTIV/ NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,IN,10,DELT,NTRANS 
COMMON /HILL/ SLOPE,COS1,COS2,SIN1,SIN2 

DIMENSION ELMAT(3,3),A(3),B8(3),XL(5),ZL(5),AL(3),BL(3),GMAT(1), 
*KZZ(1),CA(1),X(1),Z(1),H(1),NP(3,1),KRATIO(1)., 
*MATL(1),LD(1),CELMAT(3,3),BMAT(1),RHS(1),LDQ(1),QRHS(1) 

c WRITE(6,2153) 
2153 FORMAT(//,13X,’ELEM AVK (VERT.) AVCA’,//) 


Cc 
Cc 
Cc INITIALIZE GLOBAL MATRICES 
Cc 
Cc 


MUP =NNODE * (NNODE+1)/2 
DO 3 IG=1,MUP 
GMAT(1IG)=0.0 

3 BMAT(IG)=0.0 


Cc 
Cc 
C FORM DIAGONAL ADDRESSING ARRAY 
Cc 
Cc 
E=0 
DO 410 N=1,NNODE 
L=L+N 


410 LD(N)=L 
DO 300 M=1,NEL 


FORM UPPER TRIANGULAR ELEMENT MATRIX 


aan 


DO 200 1=1,3 

DO 200 J=1i,3 

ELMAT(I,J)=0.0 
200 CELMAT(I,JU)=0.0 

I=NP(1,M) 

J=NP(2,M) 


K=NP(3,M) 
Cc COMPUTE AVERAGE CONDUCTIVITY COMPONENTS FOR EACH ELEMENT 
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40 
10 


20 


30 


50 


AVKXX=0.0 

AVKZZ=0.0 

AVKXZ=0.0 

NKOUNT=1 

IF (NKOUNT-2) 10,20,30 

IA=MATL(I) 

RQ=KRATIO(IA) 

NUM=I 

GO TO 50 

IA=MATL(J) 

RQ=KRATIO(IA) 

NUM=J 

GO TO 50 

IA=MATL(K) 

RQ=KRATIO(IA) 

NUM=K 
AVKXX=AVKXX+(COS2*ROQ+SIN2) *KZZ(NUM) 
AVKZZ=AVKZZ+(SIN2*RQ+COS2) *KZZ(NUM) 
AVKXZ=AVKXZ+COS1*SIN1*KZZ(NUM) *(RQ- 1) 
AVCA=2*CA(I)+CA(JU)+CA(K) 
AVC2=CA(I1I)+2*CA(JU)+CA(K) 
AVC3=CA(I)+CA(JU)+2*CA(K) 

NKOUNT =NKOUNT + 1 

IF(NKOUNT.LE.3) GO TO 40 
AVKXX=AVKXX/3.0 

AVKZZ=AVKZZ/3.0 

AVKXZ=AVKXZ/3.0 


FORM ELEMENT DIMENSIONS 


A(1)=X(K)-X(J) 
A(2)=X(1)-X(K) 
A(3)=X(JU)-X(I) 
B(1)=Z(JU)-Z(K) 
B(2)=Z(K)-Z(1) 
B(3)=Z(1I)-Z(v) 
AREA2=A(3) *B(2)-A(2) *B(3) 


CARRY OUT INTEGRATIONS 


226 


WRITE(6, 2152) M,AVKXX,AVKZZ, AVKXZ 
FORMAT(10X,16,3E20.4) 
DEN=2.*AREA2 

WRITE(6,1000) M 


FORMAT(//,10X,2OHTHE DEN. FOR ELEMENT,13,3HIS:,//) 


WRITE(6, 1001) DEN 
FORMAT(10X,E15.6) 
Xtakatie xX 1) 
XL(2)=X(J) 
XL(3)=X(K) 
ZL(1)=Z(1) 
ZL(2)=2(u) 
ZL(3)=Z(K) 


po 245 1=1,3 
DO 246 Jv=1,3 
IF(I-J) 226,225,246 
K=(6-I-J) 
AL(1)=Xt(K)-XL(J) 
AL(2)=XL(1)-XL(K) 


256 
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BL(1)=ZL(JU)-ZL(k) 
BL(2)=ZL(K)-ZL(1) 
ELMAT(I,J)=(AVKZZ*AL(1) *AL(2)+AVKXX*BL(1)*BL(2))/DEN 
IF (SLOPE.EQ.0.0) GO TO 60 
ELMAT(1I,JU)=ELMAT(1I,J)+(AVKXZ*(AL(2) *BL(1)+BL (2) *AL(1)))/DEN 

60 CONTINUE 
GO TO 255 

225 IF (1-2) 230,235,240 

235 XL(1+2)=XL(1) 
ZL(1+2)=ZL(1) 
AVCA=AVC2 
GO TO 230 

240 XL(1+2)=XL(2) 
ZL(1+2)=ZL(2) 
AVCA=AVC3 

230 ELMAT(I,J)=(AVKZZ*((XL(1+1)-XL(1+2) ) **2)+AVKXX*((ZL(1+1)-ZL(14+2))* 
1*2))/DEN 
IF(SLOPE.EQ.0.0) GO TO 70 
ELMAT(I,JU)=ELMAT(I,JU)+(2.O0*AVKXZ*BL(1)*AL(1))/DEN 

70 CONTINUE } 

IF(NTRANS) 260,255,260 

260 CELMAT(I,J)=AVCA*AREA2/24.0 


FORM UPPER TRIANGULAR GLOBAL MATRICES,(A) AND (B) 


255 IF(I-2) 305,310,315 
305 IG=NP(1,M) 
GO TO 320 
310 IG=NP(2,M) 
GO TO 320 
315 IG=NP(3,M) 
320 IF(U-2) 325,330,335 
325 JG=NP(1,M) 
GO TO 340 
330 JG=NP(2,M) 
GO TO 340 
335 UG=NP(3,M) 
340 IF(IG-JG) 405,415,425 
425 IIG=IG 
IIJU=UG 
JG=I1G 
IG=IIJ 
405 MM=LD(UJG)-UG+IG 
GMAT (MM) =GMAT(MM)+ELMAT(I,uJ) 
GO TO 246 
415 MM=LD(JG)-JG+IG 
GMAT (MM) =GMAT(MM)+ELMAT(I,J) 
BMAT (MM) =BMAT(MM)+CELMAT(I,JU) 
246 CONTINUE 
245 CONTINUE 
WRITE(6,110) M 
WRITE(6, 250) (CELMAT(1,.J) S24, 3), F=1,3) 
WRITE(6, 250) ((CELMAT(I,JU),JU=1,3).1=1,3) 
300 CONTINUE 
WRITE(6, 345) 
WRITE(6, 354) (GMAT (MM) ,MM=1, MUP) 


WRITE(6, 355) 
WRITE(6,354) (BMAT(MM) ,MM=1,MUP) 


WIONCES wont 
)JA™ (o)uneK ae" (gj ee 


(3) dma (al leqolaew Jade ae vuaret sine 01 


of gh7n4ws3)) 
‘yam y39)) 


Cc 


a0 


COMME Qua 


aangnNAa 


OTOL 


FORM (G) MATRIX AND RHS 


2100 


650 


500 
600 


2000 
455 
a}: 
300 
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250 
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348 
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360 


IF(NTRANS) 2100,2000,2100 
DO 450 N=1,NNODE 

MM=LD(N) 

GMAT (MM) =GMAT (MM) +BMAT (MM) /DELT 
RHS(N)=0.0 


¢ 


DO 600 I=i1,NNODE 
KOLT=LD(I)-I 

KDIAG=LD(1I) 

DO 500 J=1,NNODE 

KOUNT=KOLT+J 
IF(KOUNT.GT.KDIAG) GO TO 650 
RHS(1)=RHS(1)+BMAT(KOUNT ) *H( JU) /DELT 
GO TO 500 

UM=J-1 

KROW=LD( JM) +I 
RHS(1I)=RHS(1I)+BMAT(KROW) *H(JU)/DELT 
CONTINUE 

CONTINUE 

GO TO 555 

DO 455 KS=1,NNODE 

RHS(KS)=0.0 

DO 900 L=1,MUP 

BMAT(L)=GMAT(L) 

DO 950 M=1,NNODE 

QRHS(M) =RHS(M) 

LDQ(M)=LD(M) 


PRINT OUT GLOBAL MATRICES 


WRITE(6,348) 

WRITE(6,354) (GMAT(MM) ,MM=1,MUP) 

WRITE(6,352) 

WRITE(6,354) (H(I),1=1,NNODE) 

WRITE(6,360) 

WRITE(6,354) (RHS(I),1I=1,NNODE) 

RETURN 

FORMAT(//, 10X,22HMATRIX FOR ELEMENT NO.,12,/) 

FORMAT (20X,3E25.8) 

FORMAT(//,10X,17HGLOBAL (A) MATRIX,/) 

FORMAT(//, 10X,17HGLOBAL (G) MATRIX,/) 
FORMAT(//,10X,31HTHE INITIAL SOLUTION VECTOR 1S 5/4) 
FORMAT(10E12.4) 

FORMAT(//,10X,17HGLOBAL (B) MATRIX,/) 
FORMAT(//,10X,34HTHE RHS VECTOR CT DELL~( BI Cet saz ae 


END 


SUBROUTINE HYCOND(WATCON, SLPOT, SLCON, SLCAP,KRATIO) 


THIS SUBROUTINE CALCULATES HYDRAULIC CONDUCTIVITY AND SPECIFIC 
MOISTURE CAPACITY AS A FUNCTION OF PRESSURE POTENTIAL.IT IS 
ADAPTED FROM THE WORK OF GREEN AND COREY (1971) AT THE SAVANNAH 
RIVER LABORATORY. CALCULATIONS ARE BASED ON PAPERS BY MARSHALL 


258 


Sse 0's tom ) ee a 


: Seieinasab veer 


Aes dt) op (@sraw’. FA 
f9an\40) Sumiede 


$4 tity’ 


ee Lwtray) Tasege ( iieinetidene 


uO" rates 
ane, pian 
LJ) TAMOe CS TAME 
(ny ate : a ‘ 
: tries tm y0ed me ai. 
abprwvim saegue 7 ag 


rae sae.a - 
eee: vellin, (0a) T AMD) feds 


Gabor at ct) tbe gat | 
(Oat Oi 
aginst er. Lbyendy: taRelenar 


£3), Obi, bn Anat seveons mor VA 
a (#, 290 NO 


¢ ), AMON? KOR ONY 

ARTE 2) TAD Joust. xor , VN 

NN ety naTaaw 7 y 30> ag 1 pet start. 

A fy mane ig). se 

C\\. 124 (upetayy Be + itMee, xor , \¥ 
, Ng? 


ab Tt. JATTM AROS au 
Haiavae SA ‘User a 
JJAHQaoW Va) 283984 WO_.932 


839 
840 
841 
842 
843 
844 
845 
846 
847 
848 
849 
850 
851 
852 
853 
854 
855 
856 
857 
858 
859 
860 
861 
862 
863 
864 
865 
866 
867 
868 
869 
870 
871 
872 
873 
874 
875 
876 
877 
878 
879 
880 
881 
882 
883 
884 
885 
886 
887 
888 
889 
890 
891 
892 
893 
894 
895 
896 
897 
898 


259 


AND MILLINGTON-QUIRK. THIS IS A VARIATION OF A PROGRAM DEVELOPED BY 
DR. RAY KUNZE TO CALCULATE HYDRAULIC CONDUCTIVITY OF POROUS SOLIDS 
FROM WATER RETENTION DATA. 
SEE 

GREENFR SOE. AND.) Ce "GOREY. [1974 SOL SClLeSOCH OAM: 

PROC. 35(1):3-8. 


NAANANAANAANA 


Ce eI I RII ta Rk III a aka Rf ak oak ake ake ok ak ok ake ate ate 


C INPUT VARIABLES 


THETA=WATER CONTENT FOR EACH INPUT POINT (CM**3/CM**3) 
INP=NUMBER OF INPUT DATA POINTS (IN IS LIMITED TO 20 WITH 
PRESENT FORMAT) 

NC=NUMBER OF INCREMENTED PORE CLASSES CHOSEN FOR CALCULATING 
DATA (NC IS LIMITED TO 25) 

TMAX=MAXIMUM WATER CONTENT (CM**3/CM**3) 

SCON=EXPERIMENTALLY OBTAINED SATURATED CONDUCTIVITY (CM/UNIT TIME) 

SCONR=KRATIO:RATIO OF HORIZONTAL TO VERTICAL HYDRAULIC 

CONDUCTIVITY (PRINCIPAL DIRECTIONS OF PERMEABILITY) 

DP=DESORPTION PRESSURE (CM OF WATER) 

RESWAT=ESTIMATE OF RESIDUAL(IMMOBILE) WATER 


INTERMEDIATE CALCULATED VARIABLES 


STDINC=STANDARD WATER CONTENT INCREMENT FOR CALCULATED 
VALUES(CM* *3/CM**3) 
TINC=INCREMENTED THETA (CM**3/CM**3) 
DPI=INCREMENTED DP FOR RESPECTIVE TINC (CM OF WATER) 
CLS=SQUARED RECIPROCAL OF NUMBER OF WATER CONTENT CLASSES 
(KL) 
PCH=INTERMEDIATE SUM OF PRODUCTS OF COEFFICIENTS AND 
HEADS IN CONDUCTIVITY EQUATION 
SPCH=FINAL SUM OF PRODUCTS 
ACF=CONVERSION FACTOR THAT TAKES INTO ACCOUNT TEMPERATURE 
AND GRAVITY INFLUENCES 
=4*1440*60*(SURFACE TENSION) **2/(8*DYN. VISCOSITY* 
DENSITY*GRAVITY) 
UNITS FOR VARIABLES IN ACF 
1440=MIN/DAY 
60=SEC/MIN 
SURTEN=NEWTON/M 
VISWAT=PASCAL*SEC 
DENWAT=G/CM**3 
“GRAVITY=CM/SEC**2 


OUTPUT VARIABLES 


CCAL=CALCULATED CONDUCTIVITY (CM/UNIT TIME),CALLED ’CALCULATED K ’ 
MATF=CONDUCTIVITY (CM/UNIT TIME) MATCHED WITH SATURATION VALUE, 
CALLED ’SAT MATCH’ 
FACTOR=EXPERIMENTALLY MEASURED SATURATED CONDUCTIVITY 
DIVIDED BY CALCULATED SATURATED CONDUCTIVITY 
THETA=WATER CONTENT AT UPPER END OF 
INCREMENT (CM**3/CM**3) . 
PRESSURE=DESORPTION PRESSURE (CM OF WATER) 
IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 MATF,KRATIO 
COMMON/HOOKUP/NS ,NC, RESWAT 
DIMENSION WATCON(NS, 1),SLP 
*TINC(51),CCAL(51), THETA(S1 


NMAANMANANNANAANANNANANANNNANANAANANNANAANANANANANAANANANAAANAAAAAAN|ADANAAAA 


OT(NS,1),SLCON(NS,1),SLCAP(NS,1), 
),DP(51).SPCH(S1).DP1(S1),.MATE(S1); 
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*KRATIO(NS) 
LS=1 
READ INPUT PARAMATERS AND VARIABLES 
REAO(CS,.4114)° ST. S11. S12, S19,$14,S15.516, $17,518,519 
FORMAT(10A8) 
READ(5,1116) INP, TMAX,SCON, SCONR,RESWAT 
FORMAT(15,F10.4,2£15.5,F5.3) 
READ(5,1120) SURTEN,DENWAT,VISWAT, TEMP, GRAVTY 
FORMAT(F10.6,F10.3,F12.8,2F5.1) 
NOTE ORDER OF INPUT DATA= 
THETA(1)=LOWEST WATER CONTENT 

DP(1)=HIGHEST PRESSURE (ABSOLUTE VALUE) 
READ(S, 1117) (THETA(J) , J=1, INP) 
FORMAT(20F4.3) 
READ(5,1119) (OP(J),J=1, INP) 
FORMAT(8F 10.2) 
KRATIO(LS)=SCONR 
CALCULATE CONVERSION FACTOR 
ACF=30. *SURTEN**2/(VISWAT *DENWAT*GRAVTY) 
ACF=ACF* 1440. 
CALCULATE INCREMENT SIZE 
RNC=NC 
STDINC=( TMAX-THETA( 1) )/RNC 
INITIALIZE TINC AND DPI ARRAYS 
TINC( 1) =THETA(1) 
DPI(1)=DP(1) 
NCP 1=NC+1 
INDEX I REFERS TO INCREMENTED VARIABLES 
INDEX JU REFERS TO INPUT DATA 
CALCULATE THETA INCREMENT LIMITS 
DO 1003 I=2,NCP1 
TINC(1I)=TINC(I-1)+STDINC 
CALCULATE PRESSURE INCREMENT LIMITS 
DO 1004 J=1,INP 
JO=u 
IF(THETA(J).GE.TINC(I))GO TO 1005 
CONT INUE ; 
DPI(1I)=((TINC(I)-THETA(UO-1))/(THETA(JO) -THETA(JO-1)))*(DP( uO) 
1-DP(JO-1))+O0P(JO-1) 
CONT INUE 


ADJUST DPI TO GIVE VALUES AT MIDPOINT OF INCREMENT 


DPI(NCP1)=0.0 
DO 1178 I=1,NC 

DPI(1)=(DPI(1)+DPI(I+1))*0.5 

CONTINUE 

CALCULATE ADJUSTED NUMBER OF CLASSES (ANC) CORRESPONDING TO TOTAL 
WATER CONTENT 

ANC=(TMAX-0.0)/STDINC ; 

CALCULATE SQUARED RECIPROCAL OF ’ANC 

CLS=(1.0/ANC) **2 

CALCULATE PRODUCT OF COEFFICIENT AND ‘HEAD’ TERMS FOR 
EACH PORE CLASS 

KL=NC 

DO 1176 J=1,NC 

NL=NCP 1-J 

PCH=0.0 

DO 1175 I=JU,NC 

PCH=PCH+(2*1+1-2*J) *(4./DPI (NL) )**2 

NL=NL-14 

SPCH(KL ) =PCH 
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959 Cc CORRECT POROSITY TERM WITH COTINC FUNCTION (CORRECTION 
960 C NEEDED ONLY WHEN LIMITED THETA RANGE IS USED) 
961 COTINC=TINC(NCP1)-RESWAT 
962 C CALCLATE K FOR A GIVEN WATER CONTENT AND PRESSURE 
963 CCAL (KL ) =SPCH(KL) *ACF *COTINC**2*CLS 
964 KL=KL-1 
965 14176 CONTINUE 
966 Cc CALCULATE MATCHING FACTOR | 
967 F ACTOR=SCON/CCAL (NC) 
968 C ADJUST TINC AND DPI VALUES AT UPPER LIMIT OF INCREMENTS 
969 c FOR PLOTTINGAND CALCULATE MATCHED CONDUCTIVITY AT 
370 Cc EACH WATER CONTENT 
971 DO 1179 I=1,NC 
972 -TINC(1)=TINC(I+1) 
973 DPI(I)=(DPI(1)+DPI(I+1))*0.5 
974 1179 MATF(1)=FACTOR*CCAL(I) 
975 DPI(NC)=0.0 
976 ° C PRINT OUTPUT 
977 1400 PRINT 1090,ST,ST1,ST2,ST3,S7T4,S1T5,S1T6,S1T7,S51T8,ST9 
978 1090 FORMAT(1H1,20X, 10A8/) 
379 PRINT 1180, INP, TMAX,SCON, FACTOR, SCONR 
3880 1180 FORMAT(T1,’OINP= ’,13,715,’TMAX= ’,F5.4,1T30, ’SCON= ‘,1PE10.3,T50, 
981 +’FACTOR= ’,1PE10.3,T70, ‘SCONR= ’,1PE10.3) 
982 PRINT 1121,SURTEN, DENWAT, VISWAT,RESWAT, TEMP 
983 1121 FORMAT( 1X, ’SURTEN=’,F10.6,’ DENWAT=’,F10.2,’ VISWAT=’,F12.8,’ 
984 +RESWAT=",F5.3,2X,’ TEMP= ’,F4.1,’ C’/) 
385 “PRINT 1406,GRAVTY,ACF 
386 1406. FORMAG(”’ GRAVITY= ‘,F5.41,”’ ACF= ’,1PE10.3) 
987 PRINT 1403 
988 1403 FORMAT(’OCLASS’ ,3X, ‘PRESSURE’ ,6X, ‘THETA’ ,7X, ‘CALCULATED K’ ,6X, ‘SAT 
389 + MATCH’ /2X, (1)%, 3X, (CM WATER)“, 4X, “(BY VOL)”, 7X, 7 (CM/ TIME)”, 8x, 7 
990 +(GM/TIME)’) 
991 PRINT 1404,(1,D0PI(1),TINC(I),CCAL(I),MATF(I),I=1,NC) 
392 4404. FORMAT (40 4.1053X ORF 10. 2..6X,, OPF 5.3, 7X, 1PE10..3, 7X, 1PE10.3) 
933 PRINT 1125, (JU, THETA(J) ,DP(UJU),JU=1, INP) 
994 1125 FORMAT(‘0’,12,’INPUT DATA FOR THE ABOVE OUTPUT’//1X,T4,’J’,T10, 
995 1’THETA’ , 120, ’PRESSURE’//(1X,72,13,T710,F4.3,T19,F10.2)) 
396 C SET UP STORED DATA ARRAYS 
397 DO 1185 I=1,NC 4 
9938 WATCON(LS,I)=TINC(I) 
999 SLPOT(LS,1I)=DPI(I) 
1000 SLCON(LS,I)=MATF(1) 
1001 4185 CONTINUE 
1002 LS=LS+1 
1003 TFCUS. GTENS)JEGO TO 2 
1004 GO TO 1086 
4005 2 CONTINUE 
1006 Cc COMPUTE SPECIFIC MOISTURE CAPACITY 
1007 DO. 800 1=1,NS 
1008 DO 800 J=1,NC 
1009 TE( UU NE... 1) 160" T0900 
1010 950 SLCAP(I,JU)=0.0 
1014 GO TO 800 Re 
N G 

a a3 Feo ce ie CHATCON(T 41) WATCON(1/9))/(SLPOT(y JeSEEOL ETE. 
1014 800 CONTINUE 

T 1500 ; : P 
se 1500 ee (OCLASS’ , 3X, ‘PRESSURE’ ,6X, ‘CAPACITY LOX CI) A ea 
1017 */(CM WATER)’,5X,’(1/CM) ‘) 


1018 DO 2000 K=1,NS 
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WRITE(6,2050) K 

FORMAT (//,3X,’SOIL NO.’,12,/) 

DO 2000 L=1,NC 

WRITE (6,2075) L,SLPOT(K,L),SLCAP(K,L) 
FORMAT(2X,13,3X,E10.3,4X,E10.3) 
CONTINUE 

RETURN 

END 


SUBROUTINE INPUT1 


THIS SUBROUTINE READS CONTROL DATA OF MASTER SYSTEM 


HK KK 


2 eK oo ok oe oo oe ke oie ok eo oie ke Ko oo ok ok ok oko a ok ok ie ake ak ake ok oe oo ok ok ok ok ok ok oe ok ok 


IMPLICIT REAL*8(A-H,0-Z) 


COMMON /MASTIV/ NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,IN,IO,DELT,NTRANS 
COMMON/HOOKUP/NS , NC, RESWAT 


DIMENSION HED(20) 
READ(IN, 1000) HED 
WRITE(1I0,2000) HED 


READ(IN, 1100) NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,NS,NC,DELT,NTRANS 
WRITE(I0,2100) NSUB,MXND,NUMNST,MASTR,MBEL,IDRY,NS,NC,DELT,NTRANS 


RETURN 


FORMAT STATEMENTS 


1000 
1100 
2000 
2100 


FORMAT (20A4) 

FORMAT(814,F10.5,14) 

FORMAT (20X,20A4///) 

FORMAT(’MASTER SYSTEM CONTROL VARIABLES’//, 
*’NUMBER OF SUBSYSTEMS =i 
*’/ MAXIMUM NUMBER OF INTER-BOUNDARY NODES’/ 
*’NODES IN ANY SUBSYSTEM , 
*’NUMBER OF NODES IN MASTER SYSTEM , 
*’/FLAG FOR MASTER MATRIX A Fall, 
*’ NUMBER OF GLOBAL EXTERNAL BOUNDARY CONDITIONS , 


l 
i 
Bey 

~ 
Sy 


*/DRY RUN FLAG Nath 1 he 
*’NUMBER OF MATERIALS =’ 14//, 
*’NUMBER OF PORE CLASSES =/,14//, 
*’TIME STEP SIZE =',F10.5//, 
*’PROBLEM TYPE (STEADY=0, TRANSIENT=1) =’ 14//) 


END 


SUBROUTINE INPUT3(NSTEP, IMAX,E) 
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Cc 

Cc THIS SUBROUTINE READS CONTROL DATA FOR EACH SUBSYSTEM 
Cc 

Cc 


IMPLICIT REAL*8(A-H,0-Z) 
G . 
COMMON /MASTIV/ NSUB ,MXND,NUMNST,MASTR,MBEL,IDRY,IN, 10,DELT,NTRANS 
COMMON /HILL/ SLOPE,COS1,COS2,SIN1,SIN2 


COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL ,.NFIBN,NIBNS,NSEEP,NTEMP 


READ(IN, 1000) IS,NNODE,NEL,NSTEP,IMAX,E,SLOPE,NFIBN, 
*NIBNS,NSEEP 

WRITE(IO0,2000) IS,NNODE,NEL,NSTEP,IMAX,E,SLOPE,NFIBN, - 
*NIBNS,NSEEP 

PI=3.141592653589793D0 

ANGLE=PI*SLOPE/ 180.0 

COS 1=DCOS(ANGLE) 

COS2=COS1**2 

SIN1=DSIN( ANGLE ) 

SIN2=SINi**2 

RETURN 


FORMAT STATEMENTS 
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1000 FORMAT(514,2E£12.5,314) 

2000 FORMAT(’1’,T30, ‘SUBSYSTEM NUMBER ’ ,I5 
*//,’NUMBER OF NODES : 
*//,’NUMBER OF ELEMENTS 2% 15 
*//,’NUMBER OF TIME STEPS : 
*//,’MAXIMUN NUMBER OF ITERATIONS PER TIME STEP , 
*//,’SPECIFIED DIFFERENCE CRITERION FOR ITERATIONS 
*//,’HILLSIDE SLOPE (DEGREES) 

*//,’LOCAL NUMBER OF FIRST INTER-BOUNDARY NODE 
*//,’NUMBER OF INTER-BOUNDARY NODES 
*//,’NUMBER OF SEEPAGE FACE NODES SES Je) 


Hou on 


SUBROUTINE INPUT4(X,Z,H,PSI.NP,MATL) 


Cc 

Cc THIS SEGMENT READS AND WRITES THE NODAL AND ELEMENT DATA, 
Cc TOGETHER WITH INITIAL PRESSURE POTENTIAL 

Cc 
Cc 


Se Roe SE ARK Co EO ek Ra ko EE SER Re ae ACSI oe i IS ee 
IMPLICIT REAL*8(A-H,0-Z) 
COMMON /MASTIV/ NSUB,MXND , NUMNST,MASTR,MBEL, IORY,IN,1I0,DELT,NTRANS 
: brian /SUBIV/ IS,NNODE ,NEL, ITTER, NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 


DIMENSION XC1) 201) CAD PST 1)5 
*NP(3,141),MATL(1) 

Cc READ PRELIMINARY INFORMATION 
WRITE(IO, 101) 

WRITE(IO, 115) 
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L=4 
READ(IN, 104) N,X(N),2Z(N) 
GO TO 4 

READ(IN, 104) N,X(N),2(N) 
DN=N-L 

DxX=(X(N)-X(L))/DN 
DZ=(Z(N)-Z(L))/DN 

L=L+1 

IF(N-L) 10,4,50 
X(L)=x(L-1)+DxX 
Z(L)=Z(L-1)+22 

MATL (L) =MATL(L-1) 

GO TO 15 

WRITE(10, 104) N,X(N),Z(N) 
IF(NNODE-N) 10,20,30 
WRITE(IO,105) N 

CALL EXIT 


READ(IN, 120) (H(N),N=1,NNODE) 
READ(IN,130) (MATL(N),N=1,NNODE) 


DO 150 N=1,NNODE 
PSI(N)=H(N) 
WRITE(C(IO, 106) 
WRITE(IO, 102) 


WRITE(IO,180) (N,X(N),Z(N),H(N) ,MATL(N) ,N=1,NNODE) 


C READ AND WRITE ELEMENT DATA 


C 


Cc 
Cc 
Cc 


39 


60 


BS 


56 


25 


40 


45 


WRITE(1I0, 107) 
WRITE(IO, 108) 

ML =O 

IF(ML.GE.NEL) GO TO 40 


READ(IN, 113) M,NP(1,M),NP(2,M),NP(3,M) 
WRITE(IO,113) M,NP(1,M),NP(2,M),NP(3,M) 


MM=ML + 1 
IF(MM.EQ.M) GO TO 25 
ML 1=ML+1 

IF(ML1.EQ.M) GO TO 25 
ML2=ML+2 

MLM1=ML - 1 
IF(MLM1.LE.0) GO TO 45 
DO 55 I=1,3 
NP(I,ML1)=NP(I,MLM1)+14 
IF(ML2.EQ.M) GO TO 25 
DO 56 I=1,3 
NP(I,ML2)=NP(I,ML)+1 
ML =ML2 

GO TO 60 

ML=M 

GO TO 39 

CONTINUE 

WRITE(10, 110) 


WRITE(10,113) (M,(NP(JU,M),J=1,3),M=1,NEL) 


WRITE(IO, 111) 
RETURN 


FORMAT STATEMENTS 


33 
100 


401 FORMAT(’1’,5X, ‘OUTPUT OF IN 


FORMAT (2F 10.4) 
FORMAT(313) 


PUT NODAL DATA’) 
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Cc 
Cc 


Cc 
C 
Cc 


102 FORMAT(///,10X,19H NODAL POINT OUTPUT , ///53H 
152, COORD 
FORMAT(1I6,2F12.2) 
FORMAT(1HO,28HERROR IN NODAL DATA,NODE=,14) 
FORMAT(’1’,5X,’OUTPUT OF COMPLETE NODAL DATA 
FORMAT(’1’,5X,’OUTPUT OF INPUT ELEMENT DATA 
FORMAT(///,10X,13H ELEMENT DATA///, 
124H ELEM 
FORMAT(’1’,5X, ‘OUTPUT OF COMPLETE ELEMENT DATA’,//) 
FORMAT( ‘’MLM1 IS LESS THAN OR EQUAL TO ZERO 
FORMAT(F8.2,4F8.4,313) 


104 
105 
106 
107 
108 


110 
111 
420 
113 
walks 


130 


180 FORMAT(I6,2F12.2,F14.2,4X,13) 


A SIMPLE MANAGER WHICH WORKS WITH 5 FIXED LENGTH COMMON BLOCKS, 


eZ COORD- 
120 FORMAT(10F8.2) 


FORMAT(416) 


FORMAT(///,10X,19H NODAL POINT OUTPUT,///’ 
MATL’ ,///) 


FORMAT (4012) 


END 


FUNCTION ISPAC(NAME,LENGTH,K) 


*) 
3 


®) 


NODE 


A 5-COLUMN NAME DIRECTORY AND POINTER DIRECTORY. 


2K KK eo Kt kK KK KK Kk OK KOK KK KK ie ook ok ok eK fe ok eo ok kK KK ok OK KOK 


REAL*8 NAMES, NAME 


CHECK IF NAME ALREADY EXISTS. 


ENTER NEW NAME IN DIRECTORY. 


10 
20 


30 


40 


60 


70 


COMMON /DIMCOM/ LAST1,LAST2,LAST3,LAST4,LASTS,MAXDIM, 
*NAMES(5,20),1PT(5,21),ICOM(5) 


I SPACE =LOCOM( NAME , kK) 
IF(ISPACE.EQ.0) GO TO 10 


GO TO 100 


GO TO (20,30,40,50,60) ,K 


EAST 1=UAST Tt 
LAST=LASTi1 

GOr 070 
LAST2=LAST2+1 
LAST=LAST2 
GO TO 70 
LAST3=LAST3+1 
LAST=LAST3 

GO TO 70 
LAST4=LAST4+1 
LAST=LAST4 
GO TO 70 
LASTS=LAST5+1 
LAST=LASTS 


IF(LAST.GT.MAXDIM) GO TO 200 
NAMES(K,LAST ) =NAME 
ISPACE=IPT(K,LAST) 
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IPT(K,LAST+1)=ISPACE+LENGTH 
IF((IPT(K,LAST+1)-1).GT.ICOM(K)) GO TO 300 
ISPAC=ISPACE 

RETURN 


EXITS RESULTING FROM DIAGNOSED ERRORS 


aan 


100 WRITE(6,1000) NAME 
1000 FORMAT(22H***NAME ALREADY EXISTS, 10X,A8) 
GO TO 400 
200 WRITE(6,2000) NAME,K 
2000 FORMAT(17H***TABLE OVERFLOW, 10X,A8,14) 
GO TO 400 
300 WRITE(6,3000) NAME,K,IPT(K,LAST),LENGTH 
3000 FORMAT(23H***COMMON AREA OVERFLOW ,A8,314) 
400 CALL EXIT 


G 
END 

(es 
Cc 
Cc 
Cc 

SUBROUTINE ITERAT(X,Z,H,PSI,PSIP,BMAT,RHS,QRHS,SINT,DIFF,Q,TPOT, 

*_DQ,KODE,NSPF,IMAX,E,NST,NSTEP) 
(e; 
Cc THIS SUBROUTINE CONTROLS THE ITERATION CYCLE. 
(6; 
COR OOO IIS IO IG III ICICI CII I OR ICI ICA Ik te 
( 
Cc 


IMPLICIT REAL*8(A-H,0-Z) 
COMMON /BORDER/ NEUM,NDIR,KFLAG, NBOUND(60) 
COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 


) 


DIMENSION H(1),PSI(1),PSIP(1),SINT(1),DIFF(1),KODE(1),BMAT(1), 
*QRHS(1),L0Q(1),0(1),NSPF(1),X(1),2(1), TPOT(1).RHS(1) 
IF(ITTER) 90,95,90 
95 DO 110 I=1,NNODE 
PSI(I)=0.0 
110 PSIP(I)=0.0 
J=0 
DO 600 I=1,NNODE 
IF(KODE(I).£Q.10) GO TO 600 
J=J+1 
PSIP(I)=SINT(J) 
PSI(I)=(H(1)+SINT(U))/2.0 
600 CONTINUE 
WRITE(6, 1000) 
WRITE(6, 105) (PSIP(I),1=1,NNODE) 
WRITE(6, 100) 
WRITE(6, 101) 
WRITE(6,105) (PSI(1I),I=1,NNODE) 
RITE(6, 1230) 
1230 PORMAT(// 5X, “THE RHS TERMS ARE (IT.):’.//) 
Cc WRITE(6, 105) (RHS(1I),1=1,NNODE) 
ITTER=1 
GO TO 85 
100 FORMAT(//, 
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FORMAT(10E 12.4) 

TEST=0.0 

DO 350 I=1,NNODE 
PSI(I)=0.0 

J=O 

DO 650 I=1,NNODE 
IF(KODE(I).£Q.10) GO TO 650 
J=J+14 

PSI(1I)=SINT(J) 

CONT INUE 


COMPUTE FLUXES 


140 


750 
850 


150 


2100 


CALL FLUX(X,Z,PSI,BMAT,QORHS,SINT,Q,TPOT,LDQ,KODE,NSPF) 


DO 140 I=1,NNODE 
DIFF(I)=PSI(1)-PSIP(I) 
TEST=DMAX1(TEST,DABS(DIFF(I))) 

CONT INUE 

IF(ITTER.EQ.1) GO TO 1700 
CHGE=TEST-CHEK 

IF(CHGE) 1500,2000, 2000 

IF (CHGE.GE.-O0.001) GO TO 2000 
CHEK=TEST 

WRITE(6,700) ITTER,TEST 
FORMAL(“ITTER=" , 19-5X, ’TEST=“, E42 .S) 
WRITE(6,800) 

FORMAT(//,5X, ’ DIFF PSI’,//) 
DO 750 I=1,NNODE 

WRITE(6,850) DIFF(I),PSI(1I) 
FORMAT(5X,2E£12.5) 

IF(TEST.LT.E) GO TO 50 
ITTER=ITTER+14 

IF(ITTER.GE.IMAX) GO TO 190 

DO 150 I=1,NNODE 

PSIP(I)=PSI(1) 
PSI(I)=(H(1)+PSI(1))/2.0 

GO TO 85 

CONTINUE 

WRITE(6, 120) 

WRITE(6,105) (PSI(I),1I=1,NNODE) 
WRITE(6, 3000) 

WRITE(6, 105) (Q(I),1=1,NNODE) 
ITTER=O 

DO 6O I=1,NNODE 

PSIP(I)=H(I) 

H(1I)=PSI(1I) 
PSI(1I)=PSI(1)+(PSI(1I)-PSIP(1))/2.0 
WRITE(6,500) NST 

FORMAT(/,10X, ‘END OF: TIME STEP’,15,//) 


WRITE(6, 300) : 
FORMAT(//, ‘INITIAL PRESSURE POTENTIAL - NEXT TIME STEP’ ,//) 


WRITE(6, 400) (H(I),1=1,NNODE) 


WRITE(6, 1400) ; 
FORMAT(//, ‘TOTAL HYDRAULIC POTENTIAL’ ,//) 


WRITE(6,400) (TPOT(I),1=1,NNODE) 
FORMAT(10F 12.1) 

WRITE(6, 3000) 

WRITE(6, 105) (Q(1),1=1,NNODE) 
IF(NST.NE.NSTEP) GO TO 2200 
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Cc WRITE(8,960) NNODE 
960 FORMAT(14) 

Cc DO 950 I=1,NNODE 

Cc TPOT(1I)=H(1I)+2(1) 


C 950 WRITE(8,900) I,X(I),Z(1),H(1), TPOT(I) 
900 FORMAT(14,4F12.1) 
2200 ITTER=0 
NST=NST+1 
Cc IF(NST.GT.NSTEP) GO TO 1600 
GO TO 85 
2000 WRITE(6,2050) CHGE 
GO TO 1600 
190 WRITE(6,200) 
200 FORMAT(//,’IMAX EXCEEDED’ ,//) 
1600 CALL EXIT 
85 CONTINUE 
RETURN 
1000 FORMAT(//, ‘FIRST ESTIMATE OF P(K+1):’,//) 
107- FORMAT (//; 10X, “P(K+1/2) =’ ,//) 
120 FORMAT(//,5X, ‘LATEST ESTIMATE OF P(K+1):’,//) 
2050 FORMAT(//,5X,’THE SOLUTION FOR THIS PROBLEM DOES NOT CONVERGE’ ,//, 
*5X,’THE LATEST CHANGE IS:’,F10.3,//) 
3000 FORMAT(//,5X,’THE COMPUTED FLUXES ARE’ ,//) 
END 


FUNCTION LOCOM(NAME ,K) 
LOCATES INDEX OF A GIVEN NAME IN NAMES DIRECTORY. 


SOK ROR OR ORO OO RI a a i kk oo kk kak kk ak ak a ai ae ak ar oe ie ako oie oe ako aie te ak aie a 


REAL*8 NAME ,NAMES 


COMMON /DIMCOM/ LAST1,LAST2,LAST3,LAST4,LAST5S,MAXDIM, 
*NAMES(5,20),1PT(5,21),1COM(5) 


GO TO (.10,20,30,40,50),K 
10 LAST=LAST1 
GO TO 60 
20 LAST=LAST2 
GO TO 60 
30 LAST=LAST3 
GO TO 60 
40 LAST=LAST4 
~GO TO 60 
50 LAST=LASTS 


60 IF(LAST.EQ.O) GO TO 200 
DO 100 M=1,LAST 
IF (NAMES(K,M) .NE.NAME ) GO TO 100 
LOCOM=M 
RETURN 
{OO CONTINUE 
200 LOCOM=O 


RETURN 


END 


: fi 
se : 
bind “ty <> 
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SUBROUTINE XPLOT(X,Z,H, THETA,KZZ,TPOT) 


THIS SUBROUTINE ENTERS APPROPIATE DATA INTO A FILE 
FOR SUBSEQUENT PLOTTING ON THE AEDS12 INTEGRATED GRAPHICS 
FACILITY 
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IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 KZZ 

COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 
DIMENSION X(1),2(1),H(1),THETA(1),KZZ( 1), TPOT(1) 


WRITE(8,960) NNODE 

DO 950 I=1,NNODE 

TPOT(I)=H(1I)+Z(1) 

WRITE(8,900) I,X(1I),Z(1),H(1), TPOT(I), THETA(I) ,KZZ(I).- 
FORMAT(14,4F12.1,F12.4,£15.5) 

FORMAT(14) 

CADIS CEXIT 

RETURN 

END 


SUBROUTINE SKYLIN(GMAT,LD,KOLHT) 


THIS SUBROUTINE FORMS AND STORES GLOBAL MATRIX IN ONE- 
DIMENSIONAL ARRAY, COMPUTES COLUMN LENGTHS AND ADJUSTS FOR 
SKYLINE, IT ALSO IDENTIFIES NEW LOCATIONS OF DIAGONAL ELEMENTS 


ROK IOI II RIOR II OK a II I I I aR IR I ak a IC I oka kkk a ke ak a ak ka a ak ak a ke ok ate 


IMPLICIT REAL*8(A-H,0-Z) 


COMMON /MASTIV/ NSUB,MXND,NUMNST,MASTR,MBEL,IORY,IN,IO,DELT,NTRANS 


COMMON /SUBIV/ 1S,NNODE ,NEL,ITTER,NEBEL ,NFIBN,NIBNS,NSEEP,NTEMP 
DIMENSION GMAT(1),LD(1),KOLHT(1) 


COMPUTE COLUMN LENGTHS AND ADJUST FOR SKYLINE 


320 


32:9 


SiS 


330 


KK=1 
MUM= 14 

MSUM=O 

IF (GMAT(MUM) ..£EQ.0) GO TO 315 
KOLHT (KK) =LD( KK) -MUM+ 1 
LENG=KOLHT (KK) 

KL=LD(KK) -KOLHT (KK) 

DO 325 II=1,LENG 

GMAT (MSUM+II)=GMAT(KL+II) 
MUM=LD(KK)+1 
MSUM=MSUM+KOLHT (KK) 
LD( KK) =MSUM 

KK=KK+1 

IF (NNODE-KK) 330,320,320 
MUM=MUM+ 4 

GO TO 320 

CONTINUE 

WRITE(I0,500 

PP aie a) (LD(I), 1=1,NNODE) 
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WRITE( 10,600) 

WRITE(IO,354) (KOLHT(1),I=1,NNODE) 

IALL=LD(NNODE ) 

WRITE(1I0,355) (GMAT(1I),1=1,IALL) 
355 FORMAT(10X, 10E10.4) 

RETURN 


FORMAT STATEMENTS 


SOO FORMAT(10X,25HDIAGONAL ADDRESSING ARRAY,//) 
600 FORMAT(10X,2OHCOLUMN HEIGHTS ARRAY,//) 
354 FORMAT(10X, 1015) 

END 


SUBROUTINE TABLOK(H,PSI,THETA,KZZ,CA,WATCON, SLPOT,SLCON, 
*SLCAP,MATL) 


TABLOK IS ADAPTED FROM THE SUBROUTINE OF THE SAME NAME USED IN 
THES TERM" “MODEL. 
SEE: 


HUEE. DS DeyeRee de FLUXMORE, JB] MANKIN] “ANDEC] Ee BEGOWMETCHE 
1977. TEHM: A TERRESTRIAL ECOSYSTEM HYDROLOGY MODEL. OAK RIDGE 


NATH AGe ENVIRONS "SCiewOIVa. OAK RIDGES TENNESSEE ee PUB ie 
1019. EDFB/IBP-76/8. ORNL/NSF/EATC-27. 


WATCON IS THE TABLE OF WATER CONTENTS 

SLPOT IS THE TABLE OF PRESSURE POTENTIALS 

SLCON IS THE TABLE OF CONDUCTIVITIES 

SLCAP IS THE TABLE OF SPECIFIC MOISTURE CAPACITIES 


HO KR aR aK ak I oR Rk a ek a ok ak ak a ake ok a a aR kak aR eke ade oe oe a ok 


IMPLICIT REAL*8(A-H,0-Z) 
REAL*8 KZZ 


COMMON/HOOKUP/NS,NC,RESWAT 
COMMON /SUBIV/ IS,NNODE,NEL,ITTER,NEBEL,NFIBN,NIBNS,NSEEP,NTEMP 


DIMENSION H( 1), THETA(1),PSI(1),KZZ(1),CA(1),MATL(1), 
*WATCON(NS, 1),SLPOT(NS, 1),SLCON(NS,1),SLCAP(NS, 1) 


WRITE(6, 100) ((WATCON(IL,KL),SLPOT(IL,KL),SLCON(IL,KL), 
*SLCAP(IL.KL),KL=1,NC),IL=1,NS) 
100 FORMAT //,FSr arr 10) 2, 6104952", E 10:3) 
NC ,NNODE,RES 
Me Nast Lal NCes. 14, (NNODES™.1 4. REGWAT=RE GEG 
WRITE(6,300) (PSI(IM),I1M=1,NNODE) 
300 FORMAT(//,F 14.2) 
WRITE(6,63) 
DO 6O N=1,NNODE 
J=MATL(N) 
IF(PSI(N).GE.0.0) GO TO 72 
p0 6% K=1,NC 
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1561 
1562 
1563 
1564 
1565 
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1574 
1575 
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1580 
1581 
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61 CONTINUE 
GO TO 60 

TOPPECTELTI2) GOTO) 7.4 

VALUES WITHIN TABLE RANGE : 
PROPTN=(DABS(PSI(N))-SLPOT(JU,I))/(SLPOT(JU,I-1)-SLPOT(UJ,I)) 
THETA(N) =WATCON(JU,1I)-(WATCON(U,1)-WATCON(J,I-1))*PROPTN 
KZZ(N)=SLCON(JU,1I)-(SLCON(JU,1)-SLCON(J,I-1))*PROPTN 
CA(N)=SLCAP(JU,1)-(SLCAP(JU,1)-SLCAP(JU,I-1))*PROPTN 
GO TO 60 

VALUES AT DRY END OUTSIDE TABLE RANGE 

71 PROPTN=(DABS(PSI(N))-SLPOT(J,1I))/(15000-SLPOT(UJ,I)) 
THETA(N)=WATCON(UJU,I)-(WATCON(JU,1I)-RESWAT) *PROPTN 
KZZ(N)=SLCON(JU,I)-(SLCON(U,1I) *PROPTN) 
CA(N)=0.0 
GO TO 60 

VALUES AT SATURATION 

72 THETA(N) =WATCON(J,NC) 
KZZ(N)=SLCON(UJ,NC) 
CA(N)=0.0 

60 CONTINUE 

6O WRITE(6,65) N,MATL(N),PSI(N), THETA(N) ,KZZ(N) ,CA(N) 


RETURN 

63 FORMAT(//,11X,’N MATL PSI THETA K2Z2Z’ 
= caP’,//) 

65 FORMAT(10X,213,4E£15.5) 
END 
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D. Appendix D. Example Of Input Data Set Required To Run 
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SUBFEM 
FERNOW 
PROTOTYPE SIMULATION: FERNOW EXPT. 
1 O 272 i oO 1 Py eas) 
CALVIN SILT LOAM - UPPER HORIZONS (TO 60 cm) 
3 0.455 1 RLOETO2 
73.49 ©.999 O-7OtH380 
T2oOMonS. 455 
10300 .00 339.0 0.0 
CALVIN SILT LOAM - LOWER HORIZONS (BELOW 60 cm) 
3 0.280 1 aTOE+ON 
73.49 0.399 0.011380 
-175.240.280 
7210.00 339.0 0.0 
1272 462 120 75 MSLOE-O2 
1 0.0 0.0 
8 0.0 U / LA 
9 150 60.6 
16 150. 235).6 
17 300. 12'1..2 
24 300. 296.2 
25 450 181.8 
Sr 450. 356.8 
33 600. 242.4 
40 600 417.4 
41 750 303.0 
48 Hao 478.0 
49 900. 363.6 
56 9300. 538.6 
Si 1050 424.2 
64 1050 599.2 
65 1200. 484.8 
U2 1200. 659.8 
US) 1350 545.4 
80 Ha5O) 720.4 
81 1500. 606.0 
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89 1650 666.6 
96 1650 841.6 
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113 2100. 848.5 
120 2100. 1O2G5 
121 2250 909.1 
4128 2250 1084.1 
129 2400 969.7 
ASG 2400 1144.7 
UReHTL 2550 103073 
144 2550 1iZOSe3 
145 2700 1030.9 
160 PNEISO) 1326.5 
169 3150 W272 
176 3150 1447.7 
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1508.3 

1393.9 

1568.9 

1454.5 

1629.5 

1515.1 

1690.1 

1575.7 

1750.7 

1636.3 

1811.3 

1652.1 

1827.1 

1667.8 

1842.8 

1683.6 

1858.6 

1699.4 

1874.4 

1715.1 

1890.1 

1730.9 

1905.9 
-4944.0 -4944.0 -4944.0 -4944.0 -4944.0 -4944.0 -4944.0 -4944. 
-4944.0 -4944.0 -4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180. 
-4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -4944. 
-6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -4944.0 -4944.0 -4944. 
-6180.0 -6180.0 -6180.0 -4944.0 -4944.0 -4944.0 -6180.0 -6180. 
-6180.0 -4944.0 -4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180. 
-4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -4944. 
-6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -4944.0 -4944.0 -4944. 
-6180.0 -6180.0 -6180.0 -4944.0 -4944.0 -4944.0 -6180.0 -6180. 
-6180.0 -4944.0 -4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180. 
-4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -4944. 
-6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -4944.0 -4944.0 -4944. 
-6180.0 -6180.0 -6180.0 -4944.0 -4944.0 -4944.0 -6180.0 -6180. 
-6180.0 -4944.0 -4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180. 
-4944.0 -4944.0 -6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -4944. 
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-6180.0 -7210.0 -7210.0 -7210.0 -6180.0 -6180.0 -6180.0 -6180. 
-7210.0 -7210.0 -6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -7210. 
-6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -7210.0 -7210.0 -7210. 
-6180.0 -6180.0 -6180.0 -7210.0 -7210.0 -7210.0 -6180.0 -6180. 
-6180.0 -7210.0 -7210.0 -7210.0 -6180.0 -6180.0 -6180.0 -6180. 
-7210.0 -7210.0 -6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -7210. 
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-7210.0 -7210.0 -6180.0 -6180.0 -6180.0 -6180.0 -6180.0 -7210. 
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